home *** CD-ROM | disk | FTP | other *** search
/ Magnum One / Magnum One (Mid-American Digital) (Disc Manufacturing).iso / d11 / frasrc14.arc / CALCFRAC.C < prev    next >
C/C++ Source or Header  |  1990-08-02  |  99KB  |  3,113 lines
  1. /*
  2. FRACTALS.C and CALCFRAC.C actually calculate the fractal images (well,
  3. SOMEBODY had to do it!).  The modules are set up so that all logic that
  4. is independent of any fractal-specific code is in CALCFRAC.C, and the
  5. code that IS fractal-specific is in FRACTALS.C. Original author Tim Wegner,
  6. but just about ALL the authors have contributed SOME code to this routine
  7. at one time or another, or contributed to one of the many massive
  8. restructurings.
  9.  
  10.    -------------------------------------------------------------------- */
  11.  
  12. #include <stdio.h>
  13. #include <string.h>
  14. #include <stdlib.h>
  15. #include <float.h>
  16. #include <math.h>
  17. #include <dos.h>
  18. #include <limits.h>
  19. #include <time.h>
  20. #include <stdarg.h>
  21. #include "fractint.h"
  22. #include "mpmath.h"
  23. #include "targa_lc.h"
  24.  
  25. extern struct complex initorbit;
  26. extern char useinitorbit;
  27. struct lcomplex linitorbit;
  28.  
  29. extern unsigned int decoderline[];
  30. extern unsigned int boxx[];
  31. extern int overflow;
  32. extern int bailout;
  33. extern char far plasmamessage[];
  34. long lmagnitud, llimit, llimit2, lclosenuff,l16triglim;
  35. struct complex init,tmp,old,new,saved;
  36. extern struct lcomplex lold;
  37. extern double tempsqrx,tempsqry;
  38. extern long ltempsqrx,ltempsqry;
  39. extern int biomorph;
  40. extern struct lcomplex linit;
  41. extern int basin;
  42. extern int cpu;
  43. extern int resave_flag;
  44. extern int dotmode;
  45.  
  46. int color, oldcolor, oldmax, oldmax10, row, col, passes, passnum;
  47. int iterations, invert;
  48. double f_radius,f_xcenter, f_ycenter; /* for inversion */
  49. extern double far *dx0, far *dy0;
  50. extern double far *dx1, far *dy1;
  51. long XXOne, FgOne, FgTwo, LowerLimit;
  52.  
  53. extern int LogFlag;
  54. extern unsigned char far *LogTable;
  55.  
  56. void (*plot)() = putcolor;
  57. extern int bound_trace_main();
  58.  
  59. extern double inversion[];          /* inversion radius, f_xcenter, f_ycenter */
  60. extern int      xdots, ydots;       /* coordinates of dots on the screen  */
  61. extern int      colors;             /* maximum colors available */
  62. extern int      andcolor;           /* colors-1                 */
  63. extern int      inside;             /* "inside" color to use    */
  64. extern int      outside;            /* "outside" color to use   */
  65. extern int      finattract;
  66. double          min_orbit;          /* orbit value closest to origin */
  67. int             min_index;          /* iteration of min_orbit */
  68. extern int      maxit;              /* try this many iterations */
  69. extern int      fractype;           /* fractal type */
  70. extern int      numpasses;          /* 0 = 1 pass, 1 = double pass */
  71. extern int      solidguessing;      /* 1 if solid-guessing */
  72. extern int      debugflag;          /* for debugging purposes */
  73. extern int      timerflag;          /* ditto */
  74. extern  int     diskvideo;          /* for disk-video klooges */
  75. extern int      calc_status;        /* status of calculations */
  76. extern long     calctime;           /* total calc time for image */
  77.  
  78. extern int rflag, rseed;
  79. extern int decomp[];
  80. extern int distest;
  81.  
  82. extern double   param[];            /* parameters */
  83. extern double   potparam[];         /* potential parameters */
  84. extern long     far *lx0, far *ly0; /* X, Y points */
  85. extern long     far *lx1, far *ly1; /* X, Y points */
  86. extern long     fudge;              /* fudge factor (2**n) */
  87. extern int      bitshift;           /* bit shift for fudge */
  88. extern long     delmin;             /* for periodicity checking */
  89. extern char potfile[];              /* potential filename */
  90.  
  91.  
  92. #ifndef sqr
  93. #define sqr(x) ((x)*(x))
  94. #endif
  95.  
  96. extern    int extraseg;
  97.  
  98. /* These are local but I don't want to pass them as parameters */
  99. static    unsigned char top;   /* flag to indicate top of calcfract */
  100. static    unsigned char c;
  101. static    int i;
  102.  
  103. extern double xxmin,xxmax,yymin,yymax,xx3rd,yy3rd; /* corners */
  104. extern long   xmin, xmax, ymin, ymax;              /* integer equivs */
  105. extern double plotmx1,plotmx2,plotmy1,plotmy2; /* real->screen conversion */
  106. extern double delxx,delxx2,delyy,delyy2;
  107. struct complex lambda;
  108. double deltaX, deltaY;
  109. double magnitude, rqlim, rqlim2;
  110. extern struct complex parm,parm2;
  111. int (*calctype)();
  112. double closenuff;
  113. int pixelpi; /* value of pi in pixels */
  114. FILE *fp_pot;
  115. int potflag; /* tells if continuous potential on  */
  116. unsigned long lm;               /* magnitude limit (CALCMAND) */
  117.  
  118. /* ORBIT variables */
  119. int     show_orbit;                     /* flag to turn on and off */
  120. int     orbit_ptr;                      /* pointer into save_orbit array */
  121. static int  far *save_orbit;            /* array to save orbit values */
  122. static int orbit_color=15;              /* XOR color */
  123.  
  124. int     ixstart, ixstop, iystart, iystop;       /* (for CALCMAND) start, stop here */
  125. int     symmetry;        /* symmetry flag */                                    
  126. int     calcmode;        /* for StandardFractal/calcmand:      */               
  127.                          /* 0 main rtn, 1/2 pass, ixstart etc  */               
  128.                          /* 1 do initialization only           */               
  129.                          /* 2 1pixel, reset periodicity checks */               
  130.                          /* 3 1pixel, no reset                 */               
  131. int     kbdcount;        /* avoids checking keyboard too often */
  132.  
  133. extern  int     integerfractal;         /* TRUE if fractal uses integer math */
  134.  
  135. char far *resume_info = NULL;           /* pointer to resume info if allocated */
  136. int resume_len;                         /* length of resume info */
  137. static int resume_offset;               /* offset in resume info gets */
  138. int resuming;                           /* nonzero if resuming after interrupt */
  139. int num_worklist;                       /* resume worklist for standard engine */
  140. struct workliststuff worklist[MAXCALCWORK];
  141. int xxstart,xxstop;                     /* these are same as worklist, */
  142. int yystart,yystop,yybegin;             /* declared as separate items  */
  143. int workpass,worksym;                   /* for the sake of calcmand    */
  144.  
  145. long timer_start,timer_interval;        /* timer(...) start & total */
  146.  
  147. struct complex far *dem_orbit;          /* temp used with distance estimator */
  148. static int dem_iter;                    /* number of entries in dem_orbit */
  149. double dem_delta, dem_width;
  150. #define DEMBAILOUT 100000.0
  151. #define DEMOVERFLOW 100000000000000.0
  152.  
  153. int  timer(int (*fractal)(), int argc,...);
  154. static int  calcmandorstd(int (*func)());
  155. static int  StandardCalc();
  156. extern int  calcmand();
  157. int  StandardPixel();
  158. extern int  calcmandasm();
  159. static void perform_worklist();
  160. static void setsymmetry(int,int);
  161. int  add_worklist(int,int,int,int,int,int,int);
  162. static int  combine_worklist();
  163. void tidy_worklist();
  164.  
  165. /* static vars for solidguess & its subroutines */
  166. static int maxblock,halfblock;
  167. static int guessplot;                   /* paint 1st pass row at a time?   */
  168. static int right_guess,bottom_guess;
  169. #define maxyblk 7    /* maxxblk*maxyblk*2 <= 4096, the size of "prefix" */
  170. #define maxxblk 202  /* each maxnblk is oversize by 2 for a "border" */
  171. /* next has a skip bit for each maxblock unit;
  172.    1st pass sets bit  [1]... off only if block's contents guessed;
  173.    at end of 1st pass [0]... bits are set if any surrounding block not guessed;
  174.    bits are numbered [..][y/16+1][x+1]&(1<<(y&15)) */
  175. extern unsigned int prefix[2][maxyblk][maxxblk]; /* common temp */
  176. /* size of next puts a limit of 2048 pixels across on solid guessing logic */
  177. extern char dstack[4096];               /* common temp, two put_line calls */
  178.  
  179. #define N_ATTR 8                    /* max number of attractors     */  
  180. int     attractors;                 /* number of finite attractors  */
  181.  
  182. struct complex  attr[N_ATTR];       /* finite attractor vals (f.p)  */
  183.  
  184. struct lcomplex lattr[N_ATTR];      /* finite attractor vals (int)  */
  185.  
  186.  
  187. #ifndef lsqr                                                            
  188. #define lsqr(x) (multiply((x),(x),bitshift))                            
  189. #endif                                                                  
  190.  
  191. /* -------------------------------------------------------------------- */
  192. /*              These variables are external for speed's sake only      */
  193. /* -------------------------------------------------------------------- */
  194.  
  195. extern struct lcomplex lold,lnew,lparm,lparm2;   /* added "lold" */
  196.  
  197. int periodicitycheck = 1;
  198.  
  199. extern double floatmin, floatmax;
  200.  
  201. extern int StandardFractal();
  202. int boundarytraceflag = 0;
  203.  
  204. extern int display3d;
  205.  
  206. calcfract()
  207. {
  208.    int oldnumpasses, oldsolidguessing, oldperiodicitycheck, olddistest;
  209.    int i;
  210.    oldnumpasses        = numpasses;          /* save these */
  211.    oldsolidguessing    = solidguessing;
  212.    oldperiodicitycheck = periodicitycheck;
  213.    olddistest          = distest;
  214.  
  215.    attractors = 0;          /* default to no known finite attractors  */
  216.    
  217.    display3d = 0;
  218.  
  219.    orbit_color = 15;
  220.    top     = 1;
  221.    init_misc();  /* set up some variables in parser.c */
  222.    if(fp_pot)
  223.    {
  224.       fclose(fp_pot);
  225.       fp_pot = NULL;
  226.    }
  227.  
  228.    if (invert)
  229.    {
  230.       floatmin = FLT_MIN;
  231.       floatmax = FLT_MAX;
  232.    }
  233.    calcmode = 0; /* default mode for StandardFractal, not btm/ssg */
  234.  
  235.    basin = 0;
  236.  
  237.    /* following delta values useful only for types with rotation disabled */
  238.    /* currently used only by bifurcation */
  239.    if (integerfractal)
  240.    {
  241.       distest = 0;
  242.       deltaX = (double)lx0[      1] / fudge - xxmin;
  243.       deltaY = yymax - (double)ly0[      1] / fudge;
  244.    }
  245.    else
  246.    {
  247.       deltaX = dx0[      1] - xxmin;
  248.       deltaY = yymax - dy0[      1];
  249.    }
  250.  
  251.    parm.x   = param[0];
  252.    parm.y   = param[1];
  253.    parm2.x  = param[2];
  254.    parm2.y  = param[3];
  255.  
  256.    if(fabs(potparam[0]) > 0.0)
  257.       potflag = 1;
  258.    else
  259.       potflag = 0;
  260.  
  261.    if (LogFlag)
  262.       /* ChkLogFlag(); */
  263.       SetupLogTable();
  264.  
  265.    lm = 4;                              /* CALCMAND magnitude limit */
  266.    lm = lm << bitshift;
  267.    /*
  268.          Continuous potential override (unused at the moment)
  269.             if (potparam[2] > 4 && potparam[2] < 64)
  270.                lm = potparam[2]*fudge;
  271.    */
  272.    /* set modulus bailout value if 3rd potparam set */
  273.  
  274.    if (fabs(potparam[2]) > 0.0)
  275.       rqlim = potparam[2];
  276.    else if (decomp[0] > 0 && decomp[1] > 0)
  277.       rqlim = (double)decomp[1];
  278.    else if(bailout)    /* user input bailout */
  279.       rqlim = bailout;
  280.    else if(biomorph != -1)   /* biomorph benefits from larger bailout */
  281.       rqlim = 100;
  282.    else
  283.       rqlim = fractalspecific[fractype].orbit_bailout;
  284.  
  285.    /* ORBIT stuff */
  286.    save_orbit = (int far *)((double huge *)dx0 + 4*MAXPIXELS);
  287.    show_orbit = 0;
  288.    orbit_ptr = 0;
  289.    if(colors < 16)
  290.       orbit_color = 1;
  291.  
  292.    if (integerfractal)          /* the bailout limit can't be too high here */
  293.       if (rqlim > 127.0) rqlim = 127.0;
  294.    if(inversion[0] != 0.0)
  295.    {
  296.       f_radius    = inversion[0];
  297.       f_xcenter   = inversion[1];
  298.       f_ycenter   = inversion[2];
  299.  
  300.       if (inversion[0] == AUTOINVERT)  /*  auto calc radius 1/6 screen */
  301.          inversion[0] = f_radius = min(fabs(xxmax - xxmin),
  302.              fabs(yymax - yymin)) / 6.0;
  303.  
  304.       if (invert < 2 || inversion[1] == AUTOINVERT)  /* xcenter not already set */
  305.       {
  306.          inversion[1] = f_xcenter = (xxmin + xxmax) / 2.0;
  307.          if (fabs(f_xcenter) < fabs(xxmax-xxmin) / 100)
  308.             inversion[1] = f_xcenter = 0.0;
  309.       }
  310.  
  311.       if (invert < 3 || inversion[2] == AUTOINVERT)  /* ycenter not already set */
  312.       {
  313.          inversion[2] = f_ycenter = (yymin + yymax) / 2.0;
  314.          if (fabs(f_ycenter) < fabs(yymax-yymin) / 100)
  315.             inversion[2] = f_ycenter = 0.0;
  316.       }
  317.  
  318.       invert = 3; /* so values will not be changed if we come back */
  319.    }
  320.  
  321.    if (potfile[0] != 0 || boundarytraceflag)                    /* potential file? */
  322.    {
  323.       numpasses = 0;                            /* disable dual-pass */
  324.       solidguessing = 0;                        /* disable solid-guessing */
  325.    }
  326.  
  327.    closenuff = delmin >> abs(periodicitycheck); /* for periodicity checking */
  328.    closenuff /= fudge;
  329.    rqlim2 = sqrt(rqlim);
  330.    if (integerfractal)          /* for integer routines (lambda) */
  331.    {
  332.       lparm.x = parm.x * fudge;    /* real portion of Lambda */
  333.       lparm.y = parm.y * fudge;    /* imaginary portion of Lambda */
  334.       lparm2.x = parm2.x * fudge;  /* real portion of Lambda2 */
  335.       lparm2.y = parm2.y * fudge;  /* imaginary portion of Lambda2 */
  336.       llimit = rqlim * fudge;      /* stop if magnitude exceeds this */
  337.       if (llimit <= 0) llimit = 0x7fffffff; /* klooge for integer math */
  338.       llimit2 = rqlim2 * fudge;    /* stop if magnitude exceeds this */
  339.       lclosenuff = closenuff * fudge;   /* "close enough" value */
  340.       l16triglim = 8L<<16;         /* domain limit of fast trig functions */
  341.       linitorbit.x = initorbit.x * fudge;
  342.       linitorbit.y = initorbit.y * fudge;
  343.    }
  344.    resuming = (calc_status == 2);
  345.    if (!resuming) /* free resume_info memory if any is hanging around */
  346.    {
  347.       end_resume();
  348.       resave_flag = 0;
  349.       calctime = 0;
  350.    }
  351.  
  352.    if (fractalspecific[fractype].calctype != StandardFractal
  353.        && fractalspecific[fractype].calctype != calcmand)
  354.    {
  355.       calctype = fractalspecific[fractype].calctype; /* per_image can override */
  356.       symmetry = fractalspecific[fractype].symmetry; /*   calctype & symmetry  */
  357.       plot = putcolor; /* defaults when setsymmetry not called or does nothing */
  358.       iystart = ixstart = yystart = xxstart = yybegin = 0;
  359.       iystop = yystop = ydots -1;
  360.       ixstop = xxstop = xdots -1;
  361.       calc_status = 1; /* mark as in-progress */
  362.       distest = 0; /* only standard escape time engine supports distest */
  363.       /* per_image routine is run here */
  364.       if (fractalspecific[fractype].per_image())
  365.       { /* not a stand-alone */
  366.          /* next two lines in case periodicity changed */
  367.          closenuff = delmin >> abs(periodicitycheck); /* for periodicity checking */
  368.          closenuff /= fudge;
  369.          lclosenuff = closenuff * fudge;        /* "close enough" value */
  370.          setsymmetry(symmetry,0);
  371.          timer(calctype,0); /* non-standard fractal engine */
  372.       }
  373.       if (check_key())
  374.       {
  375.          if (calc_status == 1) /* calctype didn't set this itself, */
  376.             calc_status = 3;   /* so mark it interrupted, non-resumable */
  377.       }
  378.       else
  379.          calc_status = 4; /* no key, so assume it completed */
  380.    }
  381.    else /* standard escape-time engine */
  382.       timer((int (*)())perform_worklist,0);
  383.    calctime += timer_interval;
  384.  
  385.    if (potfile[0] != 0)
  386.    {
  387.       if (calc_status == 2) /* no resume possible (at least, not yet) */
  388.          calc_status = 3;
  389.       /* potfile[0] = 0; */
  390.    }
  391.    if(LogTable)
  392.    {
  393.       farmemfree(LogTable);
  394.       LogTable = (char far *)0;
  395.    }
  396.    solidguessing    = oldsolidguessing;
  397.    numpasses        = oldnumpasses;
  398.    periodicitycheck = oldperiodicitycheck;
  399.    distest          = olddistest;
  400.    return((calc_status == 4) ? 0 : -1);
  401. }
  402.  
  403. /* Save/resume stuff:
  404.  
  405.    Engines which aren't resumable can simply ignore all this.
  406.  
  407.    Before calling the (per_image,calctype) routines (engine), calcfract sets:
  408.       "resuming" to 0 if new image, nonzero if resuming a partially done image
  409.       "calc_status" to 1
  410.    If an engine is interrupted and wants to be able to resume it must:
  411.       store whatever status info it needs to be able to resume later
  412.       set calc_status to 2 and return
  413.    If subsequently called with resuming!=0, the engine must restore status
  414.    info and continue from where it left off.
  415.  
  416.    Since the info required for resume can get rather large for some types,
  417.    it is not stored directly in save_info.  Instead, memory is dynamically
  418.    allocated as required, and stored in .fra files as a separate block.
  419.    To save info for later resume, an engine routine can use:
  420.       alloc_resume(maxsize,version)
  421.          Maxsize must be >= max bytes subsequently saved + 2; over-allocation
  422.          is harmless except for possibility of insufficient mem available;
  423.          undersize is not checked and probably causes serious misbehaviour.
  424.          Version is an arbitrary number so that subsequent revisions of the
  425.          engine can be made backward compatible.
  426.          Alloc_resume sets calc_status to 2 (resumable) if it succeeds; to 3
  427.          if it cannot allocate memory (and issues warning to user).
  428.       put_resume({bytes,&argument,} ... 0)
  429.          Can be called as often as required to store the info.
  430.          Arguments must not be far addresses.
  431.          Is not protected against calls which use more space than allocated.
  432.    To reload info when resuming, use:
  433.       start_resume()
  434.          initializes and returns version number
  435.       get_resume({bytes,&argument,} ... 0)
  436.          inverse of store_resume
  437.       end_resume()
  438.          optional, frees the memory area sooner than would happen otherwise
  439.  
  440.    Example, save info:
  441.       alloc_resume(sizeof(parmarray)+100,2);
  442.       put_resume(sizeof(int),&row, sizeof(int),&col,
  443.                  sizeof(parmarray),parmarray, 0);
  444.     restore info:
  445.       vsn = start_resume();
  446.       get_resume(sizeof(int),&row, sizeof(int),&col, 0);
  447.       if (vsn >= 2)
  448.          get_resume(sizeof(parmarray),parmarray,0);
  449.       end_resume();
  450.  
  451.    Engines which allocate a large far memory chunk of their own might
  452.    directly set resume_info, resume_len, calc_status to avoid doubling
  453.    transient memory needs by using these routines.
  454.  
  455.    StandardFractal, calcmand, solidguess, and bound_trace_main are a related
  456.    set of engines for escape-time fractals.  They use a common worklist
  457.    structure for save/resume.  Fractals using these must specify calcmand
  458.    or StandardFractal as the engine in fractalspecificinfo.
  459.    Other engines don't get btm nor ssg, don't get off-axis symmetry nor
  460.    panning (the worklist stuff), and are on their own for save/resume.
  461.  
  462.    */
  463.  
  464. int put_resume(int len, ...)
  465. {
  466.    va_list arg_marker;  /* variable arg list */
  467.    char *source_ptr;
  468.    char far *dest_ptr;
  469.    if (resume_info == NULL)
  470.       return(-1);
  471.    va_start(arg_marker,len);
  472.    dest_ptr = resume_info + resume_len;
  473.    while (len)
  474.    {
  475.       source_ptr = va_arg(arg_marker,char *);
  476.       resume_len += len;
  477.       while (--len >= 0) /* gross, but memcpy is a nogo near vs far */
  478.          *(dest_ptr++) = *(source_ptr++);
  479.       len = va_arg(arg_marker,int);
  480.    }
  481.    return(0);
  482. }
  483.  
  484. int alloc_resume(int alloclen, int version)
  485. {
  486.    if (resume_info != NULL) /* free the prior area if there is one */
  487.       farmemfree(resume_info);
  488.    if ((resume_info = farmemalloc((long)alloclen))== NULL)
  489.    {
  490.       setfortext();
  491.       printf("\n\n\n\n\nWarning - insufficient free memory to save status.");
  492.       printf("\nYou will not be able to resume calculating this image.");
  493.       printf("\n\n\nHit any key to continue.\n");
  494.       while (keypressed()) getakey(); /* flush any keyahead */
  495.       getakey();
  496.       setforgraphics();
  497.       calc_status = 3;
  498.       return(-1);
  499.    }
  500.    resume_len = 0;
  501.    put_resume(sizeof(int),&version,0);
  502.    calc_status = 2;
  503.    return(0);
  504. }
  505.  
  506. int get_resume(int len, ...)
  507. {
  508.    va_list arg_marker;  /* variable arg list */
  509.    char far *source_ptr;
  510.    char *dest_ptr;
  511.    if (resume_info == NULL)
  512.       return(-1);
  513.    va_start(arg_marker,len);
  514.    source_ptr = resume_info + resume_offset;
  515.    while (len)
  516.    {
  517.       dest_ptr = va_arg(arg_marker,char *);
  518.       resume_offset += len;
  519.       while (--len >= 0)
  520.          *(dest_ptr++) = *(source_ptr++);
  521.       len = va_arg(arg_marker,int);
  522.    }
  523.    return(0);
  524. }
  525.  
  526. int start_resume()
  527. {
  528.    int version;
  529.    if (resume_info == NULL)
  530.       return(-1);
  531.    resume_offset = 0;
  532.    get_resume(sizeof(int),&version,0);
  533.    return(version);
  534. }
  535.  
  536. end_resume()
  537. {
  538.    if (resume_info != NULL) /* free the prior area if there is one */
  539.    {
  540.       farmemfree(resume_info);
  541.       resume_info = NULL;
  542.    }
  543. }
  544.  
  545. static void perform_worklist()
  546. {
  547.    /* default setup a new worklist */
  548.    num_worklist = 1;
  549.    worklist[0].xxstart = 0;
  550.    worklist[0].yystart = worklist[0].yybegin = 0;
  551.    worklist[0].xxstop = xdots - 1;
  552.    worklist[0].yystop = ydots - 1;
  553.    worklist[0].pass = worklist[0].sym = 0;
  554.    if (resuming) /* restore worklist, if we can't the above will stay in place */
  555.    {
  556.       start_resume();
  557.       get_resume(sizeof(int),&num_worklist,sizeof(worklist),worklist,0);
  558.       end_resume();
  559.    }
  560.    while (num_worklist > 0)
  561.    {
  562.       calctype = fractalspecific[fractype].calctype; /* per_image can override */
  563.       symmetry = fractalspecific[fractype].symmetry; /*   calctype & symmetry  */
  564.       plot = putcolor; /* defaults when setsymmetry not called or does nothing */
  565.  
  566.       /* pull top entry off worklist */
  567.       ixstart = xxstart = worklist[0].xxstart;
  568.       ixstop  = xxstop  = worklist[0].xxstop;
  569.       iystart = yystart = worklist[0].yystart;
  570.       iystop  = yystop  = worklist[0].yystop;
  571.       yybegin  = worklist[0].yybegin;
  572.       workpass = worklist[0].pass;
  573.       worksym  = worklist[0].sym;
  574.       --num_worklist;
  575.       for (i=0; i<num_worklist; ++i)
  576.          worklist[i] = worklist[i+1];
  577.  
  578.       calc_status = 1; /* mark as in-progress */
  579.       fractalspecific[fractype].per_image();
  580.       closenuff = delmin >> abs(periodicitycheck); /* for periodicity checking */
  581.       closenuff /= fudge;
  582.       lclosenuff = closenuff * fudge;   /* "close enough" value */
  583.  
  584.       /* we'd have done next sooner, but per_image can set distest & rqlim */
  585.       if (distest) /* setup stuff for distance estimator */
  586.       {
  587.          rqlim = DEMBAILOUT;
  588.          dem_delta = ( sqrt( sqr(delxx) + sqr(delxx2) )  /* half a pixel width */
  589.              + sqrt( sqr(delyy) + sqr(delyy2) ) ) / 4;
  590.          dem_width = ( sqrt( sqr(xxmax-xxmin) + sqr(xx3rd-xxmin) ) * ydots/xdots
  591.              + sqrt( sqr(yymax-yymin) + sqr(yy3rd-yymin) ) ) / distest;
  592.          dem_orbit = (struct complex far *)
  593.             farmemalloc(((long)(maxit+1)) * sizeof(*dem_orbit));
  594.          if (dem_orbit == NULL)
  595.          {
  596.             setfortext();
  597.             printf("\n\n\n\nInsufficient memory.  Try reducing maximum iterations.\n");
  598.             printf("\n\nAny key to continue.\n");
  599.             buzzer(2);
  600.             getakey();
  601.             setforgraphics();
  602.             break;
  603.          }
  604.       }
  605.       setsymmetry(symmetry,1);
  606.       /* fractal routine (usually StandardFractal) is run here */
  607.       if (boundarytraceflag && (fractalspecific[fractype].flags&NOTRACE) == 0
  608.           && potfile[0] == 0)
  609.          bound_trace_main();
  610.       else if(solidguessing && (fractalspecific[fractype].flags&NOGUESS) == 0)
  611.          solidguess();
  612.       else
  613.          (*calctype)();
  614.       if (distest) /* release distance estimator work area */
  615.          farmemfree((unsigned char far *)dem_orbit);
  616.       if (check_key()) /* interrupted? */
  617.          break;
  618.    }
  619.    if (num_worklist > 0)
  620.    {  /* interrupted, resumable */
  621.       alloc_resume(sizeof(worklist)+10,1);
  622.       put_resume(sizeof(int),&num_worklist,sizeof(worklist),worklist,0);
  623.    }
  624.    else
  625.       calc_status = 4; /* completed */
  626. }
  627.  
  628. test()
  629. {
  630.    int startrow,startpass;
  631.    startrow = startpass = 0;
  632.    if (resuming)
  633.    {
  634.       start_resume();
  635.       get_resume(sizeof(int),&startrow,sizeof(int),&startpass,0);
  636.       end_resume();
  637.    }
  638.    teststart();
  639.    for (passes=startpass; passes <= numpasses ; passes++)
  640.    {
  641.       for (row = startrow; row <= iystop; row=row+1+numpasses)
  642.       {
  643.          register int col;
  644.          for (col = 0; col <= ixstop; col++)       /* look at each point on screen */
  645.          {
  646.             register color;
  647.             init.y = dy0[row]+dy1[col];
  648.             init.x = dx0[col]+dx1[row];
  649.             if(check_key())
  650.             {
  651.                testend();
  652.                alloc_resume(20,1);
  653.                put_resume(sizeof(int),&row,sizeof(int),&passes,0);
  654.                return(-1);
  655.             }
  656.             color = testpt(init.x,init.y,parm.x,parm.y,maxit,inside);
  657.             if (color >= colors) /* avoid trouble if color is 0 */
  658.                if (colors < 16)
  659.                   color &= andcolor;
  660.                else
  661.                   color = ((color-1) % andcolor) + 1; /* skip color zero */
  662.             (*plot)(col,row,color);
  663.             if(numpasses && (passes == 0))
  664.                (*plot)(col,row+1,color);
  665.          }
  666.       }
  667.       startrow = passes + 1;
  668.    }
  669.    testend();
  670.    return(0);
  671. }
  672.  
  673. /* BTM function to return the color value for a pixel.  */
  674.  
  675. static int far *LeftX  = (int far *)NULL;
  676. static int far *RightX = (int far *)NULL;
  677. static unsigned repeats;
  678.  
  679. int calc_xy(int mx, int my)
  680. {
  681.  
  682.    color = getcolor(mx,my); /* see if pixel is black */
  683.    if (color!=0)            /* pixel is NOT black so we must have already */
  684.    {                      /* calculated its color, so lets skip it      */
  685.       repeats++;            /* count successive easy ones */
  686.       return(color);
  687.    }
  688.    repeats = 0;         /* we'll have to work for this one wo reset counter */
  689.  
  690.    col = mx;
  691.    row = my;
  692.    color=(*calctype)();
  693.    return(color);
  694. } /* calc_xy function of BTM code */
  695.  
  696.  
  697. boundary_trace(int C, int R)   /* BTM main function */
  698. {
  699.    enum
  700.        {
  701.       North, East, South, West
  702.    } 
  703.    Dir;
  704.    int modeON, C_first, bcolor, low_row, iters, gcolor;
  705.    low_row = R;
  706.    modeON = 0;
  707.    Dir = East;
  708.    bcolor = color;
  709.    C_first = C;
  710.    iters = 0;
  711.    repeats = 0;
  712.  
  713.    /* main loop of BTM inside this loop the boundary is traced on screen! */
  714.    do
  715.    {
  716.       if(--kbdcount<=0)
  717.       {
  718.          if(check_key())
  719.             return(-1);
  720.          kbdcount=(cpu==386) ? 80 : 30;
  721.       }
  722.       iters++;          /* count times thru loop */
  723.  
  724.       if (C < LeftX[R])
  725.          LeftX[R]  = C; /* to aid in filling polygon later */
  726.       if (C > RightX[R])
  727.          RightX[R] = C; /* maintain left and right limits */
  728.       else
  729.          if (R==low_row)
  730.             if (C<=C_first) /* works 99.9% of time! */
  731.                break;
  732.       switch (Dir)
  733.       {
  734.       case North :
  735.          if (R > low_row)
  736.             if(calc_xy(C,R-1)==bcolor)
  737.             {
  738.                R--;
  739.                if (C > ixstart)
  740.                   if(calc_xy(C-1,R)==bcolor)
  741.                   {
  742.                      C--;
  743.                      Dir = West;
  744.                   }
  745.                break;
  746.             }
  747.          Dir = East;
  748.          break;
  749.       case East :
  750.          if (C < ixstop)
  751.             if(calc_xy(C+1,R)==bcolor)
  752.             {
  753.                C++;
  754.                if (R > low_row)
  755.                   if(calc_xy(C,R-1)==bcolor)
  756.                   {
  757.                      R--;
  758.                      Dir = North;
  759.                   }
  760.                break;
  761.             }
  762.          Dir = South;
  763.          break;
  764.       case South :
  765.          if (R < iystop)
  766.             if(calc_xy(C,R+1)==bcolor)
  767.             {
  768.                R++;
  769.                if (C < ixstop)
  770.                   if(calc_xy(C+1,R)==bcolor)
  771.                   {
  772.                      C++;
  773.                      Dir = East;
  774.                   }
  775.                break;
  776.             }
  777.          Dir = West;
  778.          break;
  779.       case West:
  780.          if (C > ixstart)
  781.             if(calc_xy(C-1,R)==bcolor)
  782.             {
  783.                C--;
  784.                if (R < iystop)
  785.                   if(calc_xy(C,R+1)== bcolor)
  786.                   {
  787.                      R++;
  788.                      Dir = South;
  789.                   }
  790.                break;
  791.             }
  792.          Dir = North;
  793.          break;
  794.       } /* case */
  795.    }
  796.    while (repeats<30000); /* emergency backstop, should never be needed */
  797.    /* PB, made above very high to allow for resumes;  did some checking
  798.          of code first, and testing, to confirm that it seems unnecessary */
  799.    if (iters<4)
  800.    {
  801.       LeftX[low_row] = 3000;
  802.       RightX[low_row] = -3000;
  803.       if (low_row+1 <= iystop)
  804.       {
  805.          LeftX[low_row+1] = 3000;
  806.          RightX[low_row+1] = -3000;
  807.       }
  808.       return(0);  /* no need to fill a polygon of 3 points! */
  809.    }
  810.  
  811.    /* Avoid tracing around whole fractal object */
  812.    if (iystop+1==ydots)
  813.       if (LeftX[0]==0)
  814.          if (RightX[0]==ixstop)
  815.             if (LeftX[iystop]==0)
  816.                if (RightX[iystop]==ixstop)
  817.                {
  818.                   /* clean up in this RARE case or next fills will fail! */
  819.                   for (low_row = 0; low_row <= ydots-1; low_row++)
  820.                   {
  821.                      LeftX[low_row] = 3000;
  822.                      RightX[low_row] = -3000;
  823.                   }
  824.                   return(0);
  825.                }
  826.    /* fill in the traced polygon, simple but it works darn well */
  827.    C = 0;
  828.    for (R = low_row; R<iystop; R++)
  829.       if (RightX[R] != -3000)
  830.       {
  831.          if((kbdcount-=2)<=0)
  832.          {
  833.             if(check_key())
  834.                return(-1);
  835.             kbdcount=(cpu==386) ? 80 : 30;
  836.          }
  837.          if(debugflag==1946)
  838.             C = fillseg1(LeftX[R], RightX[R],R, bcolor);
  839.          else
  840.             C = fillseg(LeftX[R], RightX[R],R, bcolor);
  841.  
  842.          LeftX[R]  =  3000;
  843.          RightX[R] = -3000; /* reset array element */
  844.       }
  845.       else if (C!=0) /* this is why C = 0 above! */
  846.          return(0);
  847.    return(0);
  848. } /* BTM function */
  849.  
  850. fillseg1(int LeftX, int RightX, int R,  int bcolor)
  851. {
  852.    register modeON, C;
  853.    int  gcolor;
  854.    modeON = 0;
  855.    for (C = LeftX; C <= RightX; C++)
  856.    {
  857.       gcolor=getcolor(C,R);
  858.       if (modeON!=0 && gcolor==0)
  859. /*         (*plot)(C,R,bcolor); */
  860.          (*plot)(C,R,1); /* show boundary by only filling with color 1 */
  861.       else
  862.       {
  863.          if (gcolor==bcolor) /* TW saved a getcolor here */
  864.             modeON = 1;
  865.          else
  866.             modeON = 0;
  867.       }
  868.    }
  869.    return(C);
  870. }
  871.  
  872. fillseg(int LeftX, int RightX, int R,  int bcolor)
  873. {
  874.    unsigned char *forwards;
  875.    unsigned char *backwards;
  876.    register modeON, C;
  877.    int  gcolor,i;
  878.    modeON = 0;
  879.    forwards  = (unsigned char *)decoderline;
  880.    backwards = (unsigned char *)dstack;
  881.    modeON = 0;
  882.    get_line(R,LeftX,RightX,forwards);
  883.    for (C = LeftX; C <= RightX; C++)
  884.    {
  885.       gcolor=forwards[C-LeftX];
  886.       if (modeON!=0 && gcolor==0)
  887.          forwards[C-LeftX]=bcolor;
  888.       else
  889.       {
  890.          if (gcolor==bcolor) /* TW saved a getcolor here */
  891.             modeON = 1;
  892.          else
  893.             modeON = 0;
  894.       }
  895.    }
  896.    if(plot==putcolor) /* no symmetry! easy! */
  897.       put_line(R,LeftX,RightX,forwards);
  898.    else if(plot==symplot2) /* X-axis symmetry */
  899.    {
  900.       put_line(R,   LeftX,RightX,forwards);
  901.       if ((i=yystop-(R-yystart)) > iystop)
  902.          put_line(i,LeftX,RightX,forwards);
  903.    }
  904.    else if(plot==symplot2J) /* Origin symmetry */
  905.    {
  906.       reverse_string(backwards,forwards,RightX-LeftX+1);
  907.       put_line(R,   LeftX,                  RightX,                forwards);
  908.       if ((i=yystop-(R-yystart)) > iystop)
  909.          put_line(i,xxstop-(RightX-ixstart),xxstop-(LeftX-ixstart),backwards);
  910.    }
  911.    else if(plot==symplot2Y) /* Y-axis symmetry */
  912.    {
  913.       reverse_string(backwards,forwards,RightX-LeftX+1);
  914.       put_line(R,LeftX,                  RightX,                forwards);
  915.       put_line(R,xxstop-(RightX-ixstart),xxstop-(LeftX-ixstart),backwards);
  916.    }
  917.    else if(plot==symplot4) /* X-axis and Y-axis symmetry */
  918.    {
  919.       reverse_string(backwards,forwards,RightX-LeftX+1);
  920.       put_line(R,LeftX,                     RightX,                forwards);
  921.       put_line(R,xxstop-(RightX-ixstart),   xxstop-(LeftX-ixstart),backwards);
  922.       if ((i=yystop-(R-yystart)) > iystop)
  923.       {
  924.          put_line(i,LeftX,                  RightX,                forwards);
  925.          put_line(i,xxstop-(RightX-ixstart),xxstop-(LeftX-ixstart),backwards);
  926.       }
  927.    }
  928.    else  /* the other symmetry types are on their own! */
  929.    {
  930.       int i;
  931.       for(i=LeftX;i<=RightX;i++)
  932.          (*plot)(i,R,forwards[i-LeftX]);
  933.    }
  934.    return(C);
  935. }
  936.  
  937. /* copy a string backwards for symmetry functions */
  938. reverse_string(char *t, char *s, int len)
  939. {
  940.    register i;
  941.    len--;
  942.    for(i=0;i<=len;i++)
  943.       t[i] = s[len-i];
  944.    return(0);
  945. }
  946.  
  947. bound_trace_main()
  948. {
  949.    int srow, scol;
  950.    long maxrow;
  951.    maxrow = ((long)ydots)*sizeof(int);
  952.  
  953.    if((LeftX  = (int far *)farmemalloc(maxrow))==(int far *)NULL)
  954.       return(-1);
  955.    else if((RightX  = (int far *)farmemalloc(maxrow))==(int far *)NULL)
  956.    {
  957.       farmemfree((char far *)LeftX);
  958.       return(-1);
  959.    }
  960.  
  961.    numpasses = 0; /* for calcmand */
  962.  
  963.    for (srow = 0; srow < ydots; srow++)
  964.    {
  965.       LeftX[srow] = 3000;
  966.       RightX[srow] = -3000;
  967.    }
  968.  
  969.    /* call calc routine once to get its vars initialized */
  970.    calcmode = 1;
  971.    (*calctype)();
  972.  
  973.    for (srow = iystart; srow <= iystop; srow++)
  974.    {
  975.       for (scol = ixstart; scol <= ixstop; scol++)
  976.       {
  977.          if(--kbdcount<=0)
  978.          {
  979.             if(check_key())
  980.             {
  981.                if (iystop != yystop)
  982.                   iystop = yystop - (srow - yystart); /* allow for sym */
  983.                add_worklist(xxstart,xxstop,srow,iystop,srow,0,worksym);
  984.                farmemfree((char far *)LeftX);
  985.                farmemfree((char far *)RightX);
  986.                return(-1);
  987.             }
  988.             kbdcount=(cpu==386) ? 80 : 30;
  989.          }
  990.  
  991.          /* BTM Hook! */
  992.          color = getcolor(scol,srow);
  993.          /* if pixel is BLACK (0) then we haven't done it yet!
  994.             so first calculate its color and call the routine
  995.             that will try and trace a polygon if one exists */
  996.          if (color==0)
  997.          {
  998.             calcmode = 2; /* reset periodicity checking */
  999.             row = srow;
  1000.             col = scol;
  1001.             color=(*calctype)();
  1002.             calcmode = 3; /* keep periodicity checking from pixel to pixel */
  1003.             boundary_trace(scol, srow); /* go trace boundary! WHOOOOOOO! */
  1004.          }
  1005.       }
  1006.    }
  1007.    farmemfree((char far *)LeftX);
  1008.    farmemfree((char far *)RightX);
  1009.    return(0);
  1010. } /* end of bound_trace_main */
  1011.  
  1012. StandardFractal() /* 1 or 2 pass main routine, or btm/ssg per pixel */
  1013. {
  1014.    if (calcmode > 1) /* calc 1 pixel, row/col set by caller */
  1015.    {
  1016.       if (calcmode == 2) /* reset periodicity before doing the pixel */
  1017.       {
  1018.          oldcolor = 1;
  1019.          oldmax = min(maxit, 250);
  1020.          oldmax10 = oldmax - 10;
  1021.       }
  1022.       return (StandardPixel());
  1023.    }
  1024.    kbdcount=(cpu==386) ? 80 : 30; /* one time initialization */
  1025.    if (calcmode == 1) /* init is all */
  1026.       return(0);
  1027.    return(calcmandorstd(StandardCalc));
  1028. }
  1029.  
  1030. int calcmand()
  1031. {
  1032.    if (calcmode) /* btm or ssg, 1 pixel or init call */
  1033.       return (calcmandasm());
  1034.    return (calcmandorstd(calcmandasm));
  1035. }
  1036.  
  1037. static int calcmandorstd(int (*calcrtn)())
  1038. {
  1039.    int i;
  1040.    if (numpasses && workpass == 0) /* do 1st pass of two */
  1041.    {
  1042.       passnum = 1;
  1043.       if ((*calcrtn)() == -1)
  1044.       {
  1045.          add_worklist(xxstart,xxstop,yystart,yystop,row,0,worksym);
  1046.          return(-1);
  1047.       }
  1048.       if (num_worklist > 0) /* worklist not empty, defer 2nd pass */
  1049.       {
  1050.          add_worklist(xxstart,xxstop,yystart,yystop,yystart,1,worksym);
  1051.          return(0);
  1052.       }
  1053.       workpass = 1;
  1054.       yybegin = yystart;
  1055.    }
  1056.    passnum = 2; /* second or only pass */
  1057.    if ((*calcrtn)() == -1)
  1058.    {
  1059.       i = yystop;
  1060.       if (iystop != yystop) /* must be due to symmetry */
  1061.          i -= row - iystart;
  1062.       add_worklist(xxstart,xxstop,row,i,row,workpass,worksym);
  1063.       return(-1);
  1064.    }
  1065.    return(0);
  1066. }
  1067.  
  1068. static int StandardCalc()
  1069. {
  1070.    row = yybegin;
  1071.    while (row <= iystop)
  1072.    {
  1073.       oldcolor = 1;
  1074.       oldmax = min(maxit, 250);
  1075.       oldmax10 = oldmax - 10;
  1076.       col = ixstart;
  1077.       while (col <= ixstop)
  1078.       {
  1079.          /* on 2nd pass of two, skip even pts */
  1080.          if(passnum==1 || numpasses==0 || (row&1) != 0 || (col&1) != 0)
  1081.          {
  1082.             if (StandardPixel() == -1) /* interrupted */
  1083.                return(-1);
  1084.             if (passnum == 1) /* first pass, copy pixel and bump col */
  1085.             {
  1086.                if ((row&1) == 0 && row < iystop)
  1087.                {
  1088.                   (*plot)(col,row+1,color);
  1089.                   if ((col&1) == 0 && col < ixstop)
  1090.                      (*plot)(col+1,row+1,color);
  1091.                }
  1092.                if ((col&1) == 0 && col < ixstop)
  1093.                   (*plot)(++col,row,color);
  1094.             }
  1095.          }
  1096.          ++col;
  1097.       }
  1098.       if (passnum == 1 && (row&1) == 0)
  1099.          ++row;
  1100.       ++row;
  1101.    }
  1102.    return(0);
  1103. }
  1104.  
  1105. int StandardPixel()                     /* per pixel 1/2/b/g, called with row & col set */
  1106. {
  1107.    int caught_a_cycle;
  1108.    int savedand, savedincr;     /* for periodicity checking */
  1109.    struct lcomplex lsaved;
  1110.    int i, attracted;
  1111.  
  1112.    /* really fractal specific, but we'll leave it here */
  1113.    if (!integerfractal)
  1114.    {
  1115.       saved.x = 0;
  1116.       saved.y = 0;
  1117.       init.y = dy0[row] + dy1[col];
  1118.    }
  1119.    else
  1120.    {
  1121.       lsaved.x = 0;
  1122.       lsaved.y = 0;
  1123.       linit.y = ly0[row] + ly1[col];
  1124.    }
  1125.    orbit_ptr = 0;
  1126.    color = 0;
  1127.    caught_a_cycle = 1;
  1128.    savedand = 1;                /* begin checking every other cycle */
  1129.    savedincr = 1;               /* start checking the very first time */
  1130.  
  1131.    if (inside == -61 || inside == -60)
  1132.    {
  1133.       magnitude = lmagnitud = 0;
  1134.       min_orbit = 100000.0;
  1135.    }
  1136.    overflow = 0;                /* reset integer math overflow flag */
  1137.    if (distest)
  1138.    {
  1139.       dem_iter = 0;
  1140.       if (fractalspecific[fractype].per_pixel()) /* mandels do the 1st iter */
  1141.       {
  1142.          dem_orbit[0].x = dem_orbit[0].y = 0;
  1143.          ++dem_iter;
  1144.       }
  1145.    }
  1146.    else
  1147.       fractalspecific[fractype].per_pixel(); /* initialize the calculations */
  1148.    while (++color < maxit)
  1149.    {
  1150.       attracted = FALSE;        /* check for convergence to any   */
  1151.       if (distest)
  1152.          dem_orbit[dem_iter++] = old;
  1153.       /* calculation of one orbit goes here */
  1154.       /* input in "old" -- output in "new" */
  1155.  
  1156.       if (fractalspecific[fractype].orbitcalc())
  1157.          break;
  1158.  
  1159.       if (inside == -60 || inside == -61)
  1160.       {
  1161.          if (integerfractal)
  1162.          {
  1163.             if (lmagnitud == 0)
  1164.                lmagnitud = lsqr(lnew.x) + lsqr(lnew.y);
  1165.             magnitude = lmagnitud;
  1166.             magnitude = magnitude / fudge;
  1167.          }
  1168.          else
  1169.             if (magnitude == 0.0)
  1170.                magnitude = sqr(new.x) + sqr(new.y);
  1171.          if (magnitude < min_orbit)
  1172.          {
  1173.             min_orbit = magnitude;
  1174.             min_index = color + 1;
  1175.          }
  1176.       }
  1177.       if (show_orbit)
  1178.       {
  1179.          if (!integerfractal)
  1180.             plot_orbit(new.x, new.y, -1);
  1181.          else
  1182.             iplot_orbit(lnew.x, lnew.y, -1);
  1183.       }
  1184.  
  1185.       if (attractors > 0)       /* finite attractor in the list   */
  1186.       {                         /* NOTE: Integer code is UNTESTED */
  1187.          if (integerfractal)
  1188.          {
  1189.             for (i = 0; i < attractors; i++)
  1190.             {
  1191.                if (labs(lnew.x - lattr[i].x) < lclosenuff)
  1192.                   if (labs(lnew.y - lattr[i].y) < lclosenuff)
  1193.                   {
  1194.                      attracted = TRUE;
  1195.                      break;
  1196.                   }
  1197.             }
  1198.          }
  1199.          else
  1200.          {
  1201.             for (i = 0; i < attractors; i++)
  1202.             {
  1203.                if (fabs(new.x - attr[i].x) < closenuff)
  1204.                   if (fabs(new.y - attr[i].y) < closenuff)
  1205.                   {
  1206.                      attracted = TRUE;
  1207.                      break;
  1208.                   }
  1209.             }
  1210.          }
  1211.       }
  1212.       if (attracted)
  1213.          break;                 /* AHA! Eaten by an attractor */
  1214.  
  1215.       if (oldcolor >= oldmax10)
  1216.          if (periodicitycheck)  /* only if this is OK to do */
  1217.          {
  1218.             if ((color & savedand) == 0)        /* time to save a new value */
  1219.             {
  1220.                if (!integerfractal)
  1221.                   saved = new;  /* floating pt fractals */
  1222.                else
  1223.                   lsaved = lnew;/* integer fractals */
  1224.                if (--savedincr == 0)    /* time to lengthen the periodicity? */
  1225.                {
  1226.                   savedand = (savedand << 1) + 1;       /* longer periodicity */
  1227.                   savedincr = 4;/* restart counter */
  1228.                }
  1229.             }
  1230.             else                /* check against an old save */
  1231.             {
  1232.                if (!integerfractal)     /* floating-pt periodicity chk */
  1233.                {
  1234.                   if (fabs(saved.x - new.x) < closenuff)
  1235.                      if (fabs(saved.y - new.y) < closenuff)
  1236.                      {
  1237.                         caught_a_cycle = 7;
  1238.                         color = maxit - 1;
  1239.                      }
  1240.                }
  1241.                else             /* integer periodicity check */
  1242.                {
  1243.                   if (labs(lsaved.x - lnew.x) < lclosenuff)
  1244.                      if (labs(lsaved.y - lnew.y) < lclosenuff)
  1245.                      {
  1246.                         caught_a_cycle = 7;
  1247.                         color = maxit - 1;
  1248.                      }
  1249.                }
  1250.             }
  1251.          }
  1252.    }
  1253.    if (color >= maxit)
  1254.    {
  1255.       if (oldcolor < maxit)
  1256.       {
  1257.          oldmax = oldcolor;
  1258.          oldmax10 = oldmax - 10;
  1259.       }
  1260.       if (periodicitycheck < 0 && caught_a_cycle == 7)
  1261.          color = caught_a_cycle;/* show periodicity */
  1262.    }
  1263.    if (show_orbit)
  1264.       scrub_orbit();
  1265.    oldcolor = color;
  1266.    if (color == 0)
  1267.       color = 1;                /* needed to make same as calcmand */
  1268.  
  1269.    if (distest)
  1270.    {
  1271.       double dist,temp;
  1272.       struct complex deriv;
  1273.       if (color < maxit && caught_a_cycle != 7) /* appears to be outside */
  1274.       {
  1275.          /* Distance estimator for points near Mandelbrot set */
  1276.          /* Original code by Phil Wilson, hacked around by PB */
  1277.          /* Algorithms from Peitgen & Saupe, Science of Fractal Images, p.198 */
  1278.          deriv.x = 1; /* preset and skip 1st orbit */
  1279.          deriv.y = 0;
  1280.          i = 0;
  1281.          while (++i < dem_iter)
  1282.          {
  1283.             temp = 2 * (dem_orbit[i].x * deriv.x - dem_orbit[i].y * deriv.y) + 1;
  1284.             deriv.y = 2 * (dem_orbit[i].y * deriv.x + dem_orbit[i].x * deriv.y);
  1285.             deriv.x = temp;
  1286.             if (fabs(deriv.x) > DEMOVERFLOW || fabs(deriv.y) > DEMOVERFLOW)
  1287.                break;
  1288.          }
  1289.          temp = sqr(new.x) + sqr(new.y);
  1290.          dist = log(temp) * sqrt(temp) / sqrt(sqr(deriv.x) + sqr(deriv.y));
  1291.          if (dist < dem_delta)
  1292.             color = inside;
  1293.          else if (colors == 2)
  1294.             color = !inside;
  1295.          else
  1296.             color = sqrt(dist / dem_width + 1);
  1297.       }
  1298.    }
  1299.    else if (potflag)
  1300.    {
  1301.       if (integerfractal)       /* adjust integer fractals */
  1302.       {
  1303.          new.x = lnew.x;
  1304.          new.x /= fudge;
  1305.          new.y = lnew.y;
  1306.          new.y /= fudge;
  1307.       }
  1308.       magnitude = sqr(new.x) + sqr(new.y);
  1309.  
  1310.       /* MUST call potential every pixel!! */
  1311.       color = potential(magnitude, color);
  1312.    }
  1313.    else if (decomp[0] > 0)
  1314.       decomposition();
  1315.    else if (biomorph != -1)
  1316.    {
  1317.       if (integerfractal)
  1318.       {
  1319.          if (labs(lnew.x) < llimit2 || labs(lnew.y) < llimit2)
  1320.             color = biomorph;
  1321.       }
  1322.       else
  1323.          if (fabs(new.x) < rqlim2 || fabs(new.y) < rqlim2)
  1324.             color = biomorph;
  1325.    }
  1326.    if ((kbdcount -= color) <= 0)
  1327.    {
  1328.       if (check_key())
  1329.          return (-1);
  1330.       kbdcount = (cpu == 386) ? 80 : 30;
  1331.    }
  1332.    if (oldcolor >= maxit)       /* really color, not oldcolor */
  1333.    {
  1334.       if (inside == -60)
  1335.          color = sqrt(min_orbit) * 75;
  1336.       else if (inside == -61)
  1337.          color = min_index;
  1338.       else if (inside >= 0)
  1339.          color = inside;
  1340.       else if (inside == -1)
  1341.          color = maxit;
  1342.    }
  1343.    else
  1344.       if (outside >= 0 && attracted == FALSE)           /* merge escape-time stripes to one color */
  1345.          color = outside;
  1346.  
  1347.    if (LogFlag)
  1348.       color = LogTable[color];
  1349.    if (color >= colors)         /* don't use color 0 unless is 0 from
  1350.                                        * inside/outside */
  1351.       if (colors < 16)
  1352.          color &= andcolor;
  1353.       else
  1354.          color = ((color - 1) % andcolor) + 1;  /* skip color zero */
  1355.  
  1356.    (*plot) (col, row, color);
  1357.    return (color);
  1358. }
  1359.  
  1360. decomposition()
  1361. {
  1362.    static double cos45     = 0.70710678118654750; /* cos 45     degrees */
  1363.    static double sin45     = 0.70710678118654750; /* sin 45     degrees */
  1364.    static double cos22_5   = 0.92387953251128670; /* cos 22.5   degrees */
  1365.    static double sin22_5   = 0.38268343236508980; /* sin 22.5   degrees */
  1366.    static double cos11_25  = 0.98078528040323040; /* cos 11.25  degrees */
  1367.    static double sin11_25  = 0.19509032201612820; /* sin 11.25  degrees */
  1368.    static double cos5_625  = 0.99518472667219690; /* cos 5.625  degrees */
  1369.    static double sin5_625  = 0.09801714032956060; /* sin 5.625  degrees */
  1370.    static double tan22_5   = 0.41421356237309500; /* tan 22.5   degrees */
  1371.    static double tan11_25  = 0.19891236737965800; /* tan 11.25  degrees */
  1372.    static double tan5_625  = 0.09849140335716425; /* tan 5.625  degrees */
  1373.    static double tan2_8125 = 0.04912684976946725; /* tan 2.8125 degrees */
  1374.    static double tan1_4063 = 0.02454862210892544; /* tan 1.4063 degrees */
  1375.    static long lcos45     ; /* cos 45     degrees */
  1376.    static long lsin45     ; /* sin 45     degrees */
  1377.    static long lcos22_5   ; /* cos 22.5   degrees */
  1378.    static long lsin22_5   ; /* sin 22.5   degrees */
  1379.    static long lcos11_25  ; /* cos 11.25  degrees */
  1380.    static long lsin11_25  ; /* sin 11.25  degrees */
  1381.    static long lcos5_625  ; /* cos 5.625  degrees */
  1382.    static long lsin5_625  ; /* sin 5.625  degrees */
  1383.    static long ltan22_5   ; /* tan 22.5   degrees */
  1384.    static long ltan11_25  ; /* tan 11.25  degrees */
  1385.    static long ltan5_625  ; /* tan 5.625  degrees */
  1386.    static long ltan2_8125 ; /* tan 2.8125 degrees */
  1387.    static long ltan1_4063 ; /* tan 1.4063 degrees */
  1388.    static char start=1;
  1389.    int temp = 0;
  1390.    struct lcomplex lalt;
  1391.    struct complex alt;
  1392.    if(start & integerfractal)
  1393.    {
  1394.       start = 0;
  1395.       lcos45     = cos45      *fudge;
  1396.       lsin45     = sin45      *fudge;
  1397.       lcos22_5   = cos22_5    *fudge;
  1398.       lsin22_5   = sin22_5    *fudge;
  1399.       lcos11_25  = cos11_25   *fudge;
  1400.       lsin11_25  = sin11_25   *fudge;
  1401.       lcos5_625  = cos5_625   *fudge;
  1402.       lsin5_625  = sin5_625   *fudge;
  1403.       ltan22_5   = tan22_5    *fudge;
  1404.       ltan11_25  = tan11_25   *fudge;
  1405.       ltan5_625  = tan5_625   *fudge;
  1406.       ltan2_8125 = tan2_8125  *fudge;
  1407.       ltan1_4063 = tan1_4063  *fudge;
  1408.    }
  1409.    color = 0;
  1410.    if (integerfractal) /* the only case */
  1411.    {
  1412.       if (lnew.y < 0)
  1413.       {
  1414.          temp = 2;
  1415.          lnew.y = -lnew.y;
  1416.       }
  1417.  
  1418.       if (lnew.x < 0)
  1419.       {
  1420.          ++temp;
  1421.          lnew.x = -lnew.x;
  1422.       }
  1423.  
  1424.       if (decomp[0] >= 8)
  1425.       {
  1426.          temp <<= 1;
  1427.          if (lnew.x < lnew.y)
  1428.          {
  1429.             ++temp;
  1430.             lalt.x = lnew.x; /* just */
  1431.             lnew.x = lnew.y; /* swap */
  1432.             lnew.y = lalt.x; /* them */
  1433.          }
  1434.  
  1435.          if (decomp[0] >= 16)
  1436.          {
  1437.             temp <<= 1;
  1438.             if (multiply(lnew.x,ltan22_5,bitshift) < lnew.y)
  1439.             {
  1440.                ++temp;
  1441.                lalt = lnew;
  1442.                lnew.x = multiply(lalt.x,lcos45,bitshift) +
  1443.                    multiply(lalt.y,lsin45,bitshift);
  1444.                lnew.y = multiply(lalt.x,lsin45,bitshift) -
  1445.                    multiply(lalt.y,lcos45,bitshift);
  1446.             }
  1447.  
  1448.             if (decomp[0] >= 32)
  1449.             {
  1450.                temp <<= 1;
  1451.                if (multiply(lnew.x,ltan11_25,bitshift) < lnew.y)
  1452.                {
  1453.                   ++temp;
  1454.                   lalt = lnew;
  1455.                   lnew.x = multiply(lalt.x,lcos22_5,bitshift) +
  1456.                       multiply(lalt.y,lsin22_5,bitshift);
  1457.                   lnew.y = multiply(lalt.x,lsin22_5,bitshift) -
  1458.                       multiply(lalt.y,lcos22_5,bitshift);
  1459.                }
  1460.  
  1461.                if (decomp[0] >= 64)
  1462.                {
  1463.                   temp <<= 1;
  1464.                   if (multiply(lnew.x,ltan5_625,bitshift) < lnew.y)
  1465.                   {
  1466.                      ++temp;
  1467.                      lalt = lnew;
  1468.                      lnew.x = multiply(lalt.x,lcos11_25,bitshift) +
  1469.                          multiply(lalt.y,lsin11_25,bitshift);
  1470.                      lnew.y = multiply(lalt.x,lsin11_25,bitshift) -
  1471.                          multiply(lalt.y,lcos11_25,bitshift);
  1472.                   }
  1473.  
  1474.                   if (decomp[0] >= 128)
  1475.                   {
  1476.                      temp <<= 1;
  1477.                      if (multiply(lnew.x,ltan2_8125,bitshift) < lnew.y)
  1478.                      {
  1479.                         ++temp;
  1480.                         lalt = lnew;
  1481.                         lnew.x = multiply(lalt.x,lcos5_625,bitshift) +
  1482.                             multiply(lalt.y,lsin5_625,bitshift);
  1483.                         lnew.y = multiply(lalt.x,lsin5_625,bitshift) -
  1484.                             multiply(lalt.y,lcos5_625,bitshift);
  1485.                      }
  1486.  
  1487.                      if (decomp[0] == 256)
  1488.                      {
  1489.                         temp <<= 1;
  1490.                         if (multiply(lnew.x,ltan1_4063,bitshift) < lnew.y)
  1491.                            if ((lnew.x*ltan1_4063 < lnew.y))
  1492.                               ++temp;
  1493.                      }
  1494.                   }
  1495.                }
  1496.             }
  1497.          }
  1498.       }
  1499.    }
  1500.    else /* double case */
  1501.    {
  1502.       if (new.y < 0)
  1503.       {
  1504.          temp = 2;
  1505.          new.y = -new.y;
  1506.       }
  1507.       if (new.x < 0)
  1508.       {
  1509.          ++temp;
  1510.          new.x = -new.x;
  1511.       }
  1512.       if (decomp[0] >= 8)
  1513.       {
  1514.          temp <<= 1;
  1515.          if (new.x < new.y)
  1516.          {
  1517.             ++temp;
  1518.             alt.x = new.x; /* just */
  1519.             new.x = new.y; /* swap */
  1520.             new.y = alt.x; /* them */
  1521.          }
  1522.          if (decomp[0] >= 16)
  1523.          {
  1524.             temp <<= 1;
  1525.             if (new.x*tan22_5 < new.y)
  1526.             {
  1527.                ++temp;
  1528.                alt = new;
  1529.                new.x = alt.x*cos45 + alt.y*sin45;
  1530.                new.y = alt.x*sin45 - alt.y*cos45;
  1531.             }
  1532.  
  1533.             if (decomp[0] >= 32)
  1534.             {
  1535.                temp <<= 1;
  1536.                if (new.x*tan11_25 < new.y)
  1537.                {
  1538.                   ++temp;
  1539.                   alt = new;
  1540.                   new.x = alt.x*cos22_5 + alt.y*sin22_5;
  1541.                   new.y = alt.x*sin22_5 - alt.y*cos22_5;
  1542.                }
  1543.  
  1544.                if (decomp[0] >= 64)
  1545.                {
  1546.                   temp <<= 1;
  1547.                   if (new.x*tan5_625 < new.y)
  1548.                   {
  1549.                      ++temp;
  1550.                      alt = new;
  1551.                      new.x = alt.x*cos11_25 + alt.y*sin11_25;
  1552.                      new.y = alt.x*sin11_25 - alt.y*cos11_25;
  1553.                   }
  1554.  
  1555.                   if (decomp[0] >= 128)
  1556.                   {
  1557.                      temp <<= 1;
  1558.                      if (new.x*tan2_8125 < new.y)
  1559.                      {
  1560.                         ++temp;
  1561.                         alt = new;
  1562.                         new.x = alt.x*cos5_625 + alt.y*sin5_625;
  1563.                         new.y = alt.x*sin5_625 - alt.y*cos5_625;
  1564.                      }
  1565.  
  1566.                      if (decomp[0] == 256)
  1567.                      {
  1568.                         temp <<= 1;
  1569.                         if ((new.x*tan1_4063 < new.y))
  1570.                            ++temp;
  1571.                      }
  1572.                   }
  1573.                }
  1574.             }
  1575.          }
  1576.       }
  1577.    }
  1578.    for (i = 1; temp > 0; ++i)
  1579.    {
  1580.       if (temp & 1)
  1581.          color = (1 << i) - 1 - color;
  1582.       temp >>= 1;
  1583.    }
  1584.    if (decomp[0] == 2)
  1585.       color &= 1;
  1586.    if (colors > decomp[0])
  1587.       color++;
  1588. }
  1589.  
  1590. #if 0
  1591. MainNewton()
  1592. {
  1593.    Newton(); /* Lee's Newton */
  1594.    return(0);
  1595. }
  1596. #endif
  1597.  
  1598. /* Symmetry plot for period PI */
  1599. void symPIplot(x, y, color)
  1600. int x, y, color ;
  1601. {
  1602.    while(x <= xxstop)
  1603.    {
  1604.       putcolor(x, y, color) ;
  1605.       x += pixelpi;
  1606.    }
  1607. }
  1608. /* Symmetry plot for period PI plus Origin Symmetry */
  1609. void symPIplot2J(x, y, color)
  1610. int x, y, color ;
  1611. {
  1612.    int i,j;
  1613.    while(x <= xxstop)
  1614.    {
  1615.       putcolor(x, y, color) ;
  1616.       if ((i=yystop-(y-yystart)) > iystop && i < ydots
  1617.           && (j=xxstop-(x-xxstart)) < xdots)
  1618.          putcolor(j, i, color) ;
  1619.       x += pixelpi;
  1620.    }
  1621. }
  1622. /* Symmetry plot for period PI plus Both Axis Symmetry */
  1623. void symPIplot4J(x, y, color)
  1624. int x, y, color ;
  1625. {
  1626.    int i,j;
  1627.    while(x <= (xxstart+xxstop)/2)
  1628.    {
  1629.       j = xxstop-(x-xxstart);
  1630.       putcolor(       x , y , color) ;
  1631.       if (j < xdots)
  1632.          putcolor(j , y , color) ;
  1633.       if ((i=yystop-(y-yystart)) > iystop && i < ydots)
  1634.       {
  1635.          putcolor(x , i , color) ;
  1636.          if (j < xdots)
  1637.             putcolor(j , i , color) ;
  1638.       }
  1639.       x += pixelpi;
  1640.    }
  1641. }
  1642.  
  1643. /* Symmetry plot for X Axis Symmetry */
  1644. void symplot2(x, y, color)
  1645. int x, y, color ;
  1646. {
  1647.    int i;
  1648.    putcolor(x, y, color) ;
  1649.    if ((i=yystop-(y-yystart)) > iystop && i < ydots)
  1650.       putcolor(x, i, color) ;
  1651. }
  1652.  
  1653. /* Symmetry plot for Y Axis Symmetry */
  1654. void symplot2Y(x, y, color)
  1655. int x, y, color ;
  1656. {
  1657.    int i;
  1658.    putcolor(x, y, color) ;
  1659.    if ((i=xxstop-(x-xxstart)) < xdots)
  1660.       putcolor(i, y, color) ;
  1661. }
  1662.  
  1663. /* Symmetry plot for Origin Symmetry */
  1664. void symplot2J(x, y, color)
  1665. int x, y, color ;
  1666. {
  1667.    int i,j;
  1668.    putcolor(x, y, color) ;
  1669.    if ((i=yystop-(y-yystart)) > iystop && i < ydots
  1670.        && (j=xxstop-(x-xxstart)) < xdots)
  1671.       putcolor(j, i, color) ;
  1672. }
  1673.  
  1674. /* Symmetry plot for Both Axis Symmetry */
  1675. void symplot4(x, y, color)
  1676. int x, y, color ;
  1677. {
  1678.    int i,j;
  1679.    j = xxstop-(x-xxstart);
  1680.    putcolor(       x , y, color) ;
  1681.    if (j < xdots)
  1682.       putcolor(    j , y, color) ;
  1683.    if ((i=yystop-(y-yystart)) > iystop && i < ydots)
  1684.    {
  1685.       putcolor(    x , i, color) ;
  1686.       if (j < xdots)
  1687.          putcolor(j , i, color) ;
  1688.    }
  1689. }
  1690.  
  1691. /* Symmetry plot for X Axis Symmetry - Striped Newtbasin version */
  1692. void symplot2basin(x, y, color)
  1693. int x, y, color ;
  1694. {
  1695.    int i,stripe;
  1696.    extern int degree;
  1697.    putcolor(x, y, color) ;
  1698.    if(basin==2 && color > 8)
  1699.       stripe=8;
  1700.    else
  1701.       stripe = 0;
  1702.    if ((i=yystop-(y-yystart)) > iystop && i < ydots)
  1703.    {
  1704.       color -= stripe;                    /* reconstruct unstriped color */
  1705.       color = (degree+1-color)%degree+1;  /* symmetrical color */
  1706.       color += stripe;                    /* add stripe */
  1707.       putcolor(x, i,color)  ;
  1708.    }
  1709. }
  1710.  
  1711. /* Symmetry plot for Both Axis Symmetry  - Newtbasin version */
  1712. void symplot4basin(x, y, color)
  1713. int x, y, color ;
  1714. {
  1715.    extern int degree;
  1716.    int i,j,color1,stripe;
  1717.    if(color == 0) /* assumed to be "inside" color */
  1718.    {
  1719.       symplot4(x, y, color);
  1720.       return;
  1721.    }
  1722.    if(basin==2 && color > 8)
  1723.       stripe = 8;
  1724.    else
  1725.       stripe = 0;
  1726.    color -= stripe;               /* reconstruct unstriped color */
  1727.    color1 = degree/2+degree+2 - color;
  1728.    if(color < degree/2+2)
  1729.       color1 = degree/2+2 - color;
  1730.    else
  1731.       color1 = degree/2+degree+2 - color;
  1732.    j = xxstop-(x-xxstart);
  1733.    putcolor(       x, y, color+stripe) ;
  1734.    if (j < xdots)
  1735.       putcolor(    j, y, color1+stripe) ;
  1736.    if ((i=yystop-(y-yystart)) > iystop && i < ydots)
  1737.    {
  1738.       putcolor(    x, i, stripe + (degree+1 - color)%degree+1) ;
  1739.       if (j < xdots)
  1740.          putcolor(j, i, stripe + (degree+1 - color1)%degree+1) ;
  1741.    }
  1742. }
  1743.  
  1744.  
  1745.  
  1746. /* Do nothing plot!!! */
  1747. void noplot(int x,int y,int color)
  1748. {
  1749. }
  1750.  
  1751. static int iparmx;                                  /* iparmx = parm.x * 16 */
  1752. static int shiftvalue;                          /* shift based on #colors */
  1753.  
  1754. typedef struct palett
  1755. {
  1756.    unsigned char red;
  1757.    unsigned char green;
  1758.    unsigned char blue;
  1759. } 
  1760. Palettetype;
  1761.  
  1762. extern int colors;
  1763. extern int xdots,ydots;
  1764. extern Palettetype dacbox[256];
  1765. extern int cpu, daccount, cyclelimit;
  1766. extern int resave_flag;
  1767. static int plasma_check;                        /* to limit kbd checking */
  1768.  
  1769. void adjust(int xa,int ya,int x,int y,int xb,int yb)
  1770. {
  1771.    long pseudorandom;
  1772.    if(getcolor(x,y) != 0)
  1773.       return;
  1774.  
  1775.    pseudorandom = ((long)iparmx)*((rand()-16383));
  1776.    pseudorandom = pseudorandom*(abs(xa-xb)+abs(ya-yb));
  1777.    pseudorandom = pseudorandom >> shiftvalue;
  1778.    pseudorandom = ((getcolor(xa,ya)+getcolor(xb,yb)+1)>>1)+pseudorandom;
  1779.    if (pseudorandom <   1) pseudorandom =   1;
  1780.    if (pseudorandom >= colors) pseudorandom = colors-1;
  1781.    putcolor(x,y,(int)pseudorandom);
  1782. }
  1783.  
  1784. void subDivide(int x1,int y1,int x2,int y2)
  1785. {
  1786.    int x,y;
  1787.    int v;
  1788.    if ((++plasma_check & 0x7f) == 1)
  1789.       if(check_key())
  1790.       {
  1791.          plasma_check--;
  1792.          return;
  1793.       }
  1794.    if(x2-x1<2 && y2-y1<2)
  1795.       return;
  1796.    x = (x1+x2)>>1;
  1797.    y = (y1+y2)>>1;
  1798.    adjust(x1,y1,x ,y1,x2,y1);
  1799.    adjust(x2,y1,x2,y ,x2,y2);
  1800.    adjust(x1,y2,x ,y2,x2,y2);
  1801.    adjust(x1,y1,x1,y ,x1,y2);
  1802.    if(getcolor(x,y) == 0)
  1803.    {
  1804.       v = (getcolor(x1,y1)+getcolor(x2,y1)+getcolor(x2,y2)+
  1805.           getcolor(x1,y2)+2)>>2;
  1806.       putcolor(x,y,v);
  1807.    }
  1808.    subDivide(x1,y1,x ,y);
  1809.    subDivide(x ,y1,x2,y);
  1810.    subDivide(x ,y ,x2,y2);
  1811.    subDivide(x1,y ,x ,y2);
  1812. }
  1813.  
  1814. plasma()
  1815. {
  1816.    plasma_check = 0;
  1817.  
  1818.    if(colors < 4) {
  1819.       setvideomode(3,0,0,0);
  1820.       buzzer(2);
  1821.       helpmessage(plasmamessage);
  1822.       return(-1);
  1823.    }
  1824.    iparmx = param[0] * 8;
  1825.    if (parm.x <= 0.0) iparmx = 16;
  1826.    if (parm.x >= 100) iparmx = 800;
  1827.  
  1828.    if (!rflag) rseed = (int)time(NULL);
  1829.    srand(rseed);
  1830.  
  1831.    if (colors == 256)                   /* set the (256-color) palette */
  1832.       set_Plasma_palette();             /* skip this if < 256 colors */
  1833.  
  1834.    if (colors > 16)                     /* adjust based on # of colors */
  1835.       shiftvalue = 18;
  1836.    else {
  1837.       if (colors > 4)
  1838.          shiftvalue = 22;
  1839.       else {
  1840.          if (colors > 2)
  1841.             shiftvalue = 24;
  1842.          else
  1843.             shiftvalue = 25;
  1844.       }
  1845.    }
  1846.  
  1847.    putcolor(      0,      0,1+(((rand()/colors)*(colors-1))>>(shiftvalue-11)));
  1848.    putcolor(xdots-1,      0,1+(((rand()/colors)*(colors-1))>>(shiftvalue-11)));
  1849.    putcolor(xdots-1,ydots-1,1+(((rand()/colors)*(colors-1))>>(shiftvalue-11)));
  1850.    putcolor(      0,ydots-1,1+(((rand()/colors)*(colors-1))>>(shiftvalue-11)));
  1851.  
  1852.    subDivide(0,0,xdots-1,ydots-1);
  1853.    if (! check_key())
  1854.       return(0);
  1855.    else
  1856.       return(1);
  1857. }
  1858.  
  1859. set_Plasma_palette()
  1860. {
  1861.    static Palettetype Red    = {
  1862.       63, 0, 0   };
  1863.    static Palettetype Green  = { 
  1864.       0,63, 0   };
  1865.    static Palettetype Blue   = { 
  1866.       0, 0,63   };
  1867.    int i;
  1868.  
  1869.    if(CustomLut()) return(0);             /* TARGA 3 June 89 j mclain */
  1870.  
  1871.    dacbox[0].red  = 0 ;
  1872.    dacbox[0].green= 0 ;
  1873.    dacbox[0].blue = 0 ;
  1874.    for(i=1;i<=85;i++)
  1875.    {
  1876.       dacbox[i].red   = (i*Green.red   + (86-i)*Blue.red)/85;
  1877.       dacbox[i].green = (i*Green.green + (86-i)*Blue.green)/85;
  1878.       dacbox[i].blue  = (i*Green.blue  + (86-i)*Blue.blue)/85;
  1879.  
  1880.       dacbox[i+85].red   = (i*Red.red   + (86-i)*Green.red)/85;
  1881.       dacbox[i+85].green = (i*Red.green + (86-i)*Green.green)/85;
  1882.       dacbox[i+85].blue  = (i*Red.blue  + (86-i)*Green.blue)/85;
  1883.  
  1884.       dacbox[i+170].red   = (i*Blue.red   + (86-i)*Red.red)/85;
  1885.       dacbox[i+170].green = (i*Blue.green + (86-i)*Red.green)/85;
  1886.       dacbox[i+170].blue  = (i*Blue.blue  + (86-i)*Red.blue)/85;
  1887.    }
  1888.    ValidateLuts(NULL);     /* TARGA 3 June 89  j mclain */
  1889.    if (dotmode != 11)
  1890.       spindac(0,1);
  1891. }
  1892.  
  1893. check_key()
  1894. {
  1895.    int key;
  1896.    if((key = keypressed()) != 0) {
  1897.       if(key != 'o' && key != 'O')
  1898.          return(-1);
  1899.       getakey();
  1900.       if (dotmode != 11)
  1901.          show_orbit = 1 - show_orbit;
  1902.    }
  1903.    return(0);
  1904. }
  1905.  
  1906. #define RANDOM(x)  (rand()%(x))
  1907.  
  1908. /* This constant assumes that rand() returns a value from 0 to 32676 */
  1909. #define FOURPI  (long)(4*PI*(double)(1L << 16))
  1910.  
  1911. diffusion()
  1912. {
  1913.    int xmax,ymax,xmin,ymin;     /* Current maximum coordinates */
  1914.    int border;   /* Distance between release point and fractal */
  1915.    int i;
  1916.    double cosine,sine,angle;
  1917.    long lcosine,lsine;
  1918.    register int x,y;
  1919.    extern char floatflag;
  1920.  
  1921.  
  1922.    if (diskvideo)
  1923.    {
  1924.       setvideomode(3,0,0,0);
  1925.       buzzer(2);
  1926.       helpmessage(plasmamessage);
  1927.       return(-1);
  1928.    }
  1929.  
  1930.    bitshift = 16;
  1931.    fudge = 1L << 16;
  1932.  
  1933.    border = param[0];
  1934.  
  1935.    if (border <= 0)
  1936.       border = 10;
  1937.  
  1938.    if (!rflag)
  1939.       rseed = (int)time(NULL);
  1940.    srand(rseed);
  1941.  
  1942.    xmax = xdots / 2 + border;  /* Initial box */
  1943.    xmin = xdots / 2 - border;
  1944.    ymax = ydots / 2 + border;
  1945.    ymin = ydots / 2 - border;
  1946.  
  1947.    putcolor(xdots / 2, ydots / 2,RANDOM(colors-1)+1);  /* Seed point */
  1948.  
  1949.    while (1)
  1950.    {
  1951.       /* Release new point on circle just inside the box */
  1952.  
  1953.       if (floatflag)
  1954.       {
  1955.          angle=2*(double)rand()/(RAND_MAX/PI);
  1956.          FPUsincos(&angle,&sine,&cosine);
  1957.          x = cosine*(xmax-xmin) + xdots;
  1958.          y = sine  *(ymax-ymin) + ydots;
  1959.       }
  1960.       else
  1961.       {
  1962.          SinCos086(multiply((long)rand(),FOURPI,16),&lsine,&lcosine);
  1963.          x = (lcosine*(long)(xmax-xmin) >> 16) + xdots;
  1964.          y = (lsine  *(long)(ymax-ymin) >> 16) + ydots;
  1965.       }
  1966.  
  1967.       x = x >> 1; /* divide by 2 */
  1968.       y = y >> 1;
  1969.  
  1970.       /* Loop as long as the point (x,y) is surrounded by color 0 */
  1971.       /* on all eight sides                                       */
  1972.       while((getcolor(x+1,y+1) == 0) && (getcolor(x+1,y) == 0) &&
  1973.           (getcolor(x+1,y-1) == 0) && (getcolor(x  ,y+1) == 0) &&
  1974.           (getcolor(x  ,y-1) == 0) && (getcolor(x-1,y+1) == 0) &&
  1975.           (getcolor(x-1,y) == 0) && (getcolor(x-1,y-1) == 0))
  1976.       {
  1977.          /* Erase moving point */
  1978.          if (show_orbit)
  1979.             putcolor(x,y,0);
  1980.  
  1981.          /* Make sure point is inside the box */
  1982.          if (x==xmax)
  1983.             x--;
  1984.          else if (x==xmin)
  1985.             x++;
  1986.          if (y==ymax)
  1987.             y--;
  1988.          else if (y==ymin)
  1989.             y++;
  1990.  
  1991.          /* Take one random step */
  1992.          x += RANDOM(3) - 1;
  1993.          y += RANDOM(3) - 1;
  1994.  
  1995.          /* Check keyboard */
  1996.          if ((++plasma_check & 0x7f) == 1)
  1997.             if(check_key())
  1998.             {
  1999.                plasma_check--;
  2000.                return 1;
  2001.             }
  2002.  
  2003.          /* Show the moving point */
  2004.          if (show_orbit)
  2005.             putcolor(x,y,RANDOM(colors-1)+1);
  2006.  
  2007.       }
  2008.       putcolor(x,y,RANDOM(colors-1)+1);
  2009.  
  2010.       /* Is point to close to the edge? */
  2011.       if (((x+border)>xmax) || ((x-border)<xmin)
  2012.           || ((y-border)<ymin) || ((y+border)>ymax))
  2013.       {
  2014.          /* Increase box size, but not past the edge of the screen */
  2015.          if (ymin != 1)
  2016.          {
  2017.             ymin--;
  2018.             ymax++;
  2019.          }
  2020.          if (xmin != 1)
  2021.          {
  2022.             xmin--;
  2023.             xmax++;
  2024.          }
  2025.          if ((ymin==1) || (xmin==1))
  2026.             return 0;
  2027.       }
  2028.    }
  2029. }
  2030.  
  2031. /* Showing orbit requires converting real co-ords to screen co-ords.
  2032.    Define:
  2033.        Xs == xxmax-xx3rd               Ys == yy3rd-yymax
  2034.        W  == xdots-1                   D  == ydots-1
  2035.    We know that:
  2036.        realx == lx0[col] + lx1[row]
  2037.        realy == ly0[row] + ly1[col]
  2038.        lx0[col] == (col/width) * Xs + xxmin
  2039.        lx1[row] == row * delxx
  2040.        ly0[row] == (row/D) * Ys + yymax
  2041.        ly1[col] == col * (0-delyy)
  2042.   so:
  2043.        realx == (col/W) * Xs + xxmin + row * delxx
  2044.        realy == (row/D) * Ys + yymax + col * (0-delyy)
  2045.   and therefore:
  2046.        row == (realx-xxmin - (col/W)*Xs) / Xv    (1)
  2047.        col == (realy-yymax - (row/D)*Ys) / Yv    (2)
  2048.   substitute (2) into (1) and solve for row:
  2049.        row == ((realx-xxmin)*(0-delyy2)*W*D - (realy-yymax)*Xs*D)
  2050.                       / ((0-delyy2)*W*delxx2*D-Ys*Xs)
  2051.   */
  2052.  
  2053. static void plotdorbit(double dx, double dy, int color)
  2054. {
  2055.    int i, j, c;
  2056.    if (orbit_ptr >= 1500) return;
  2057.    i = dy * plotmx1 - dx * plotmx2;
  2058.    if (i < 0 || i >= xdots) return;
  2059.    j = dx * plotmy1 - dy * plotmy2;
  2060.    if (j < 0 || j >= ydots) return;
  2061.    /* save orbit value */
  2062.    if(color == -1)
  2063.    {
  2064.       *(save_orbit + orbit_ptr++) = i;
  2065.       *(save_orbit + orbit_ptr++) = j;
  2066.       *(save_orbit + orbit_ptr++) = c = getcolor(i,j);
  2067.       putcolor(i,j,c^orbit_color);
  2068.    }
  2069.    else
  2070.       putcolor(i,j,color);
  2071. }
  2072.  
  2073. iplot_orbit(ix, iy, color)
  2074. long ix, iy;
  2075. int color;
  2076. {
  2077.    plotdorbit((double)ix/fudge-xxmin,(double)iy/fudge-yymax,color);
  2078. }
  2079.  
  2080. plot_orbit(real,imag,color)
  2081. double real,imag;
  2082. int color;
  2083. {
  2084.    plotdorbit(real-xxmin,imag-yymax,color);
  2085. }
  2086.  
  2087. scrub_orbit()
  2088. {
  2089.    int i,j,c;
  2090.    while(orbit_ptr > 0)
  2091.    {
  2092.       c = *(save_orbit + --orbit_ptr);
  2093.       j = *(save_orbit + --orbit_ptr);
  2094.       i = *(save_orbit + --orbit_ptr);
  2095.       putcolor(i,j,c);
  2096.    }
  2097. }
  2098.  
  2099. /* timer function. Assumes #include <time.h> */
  2100. int timer(int (*fractal)(),int argc,...)
  2101. {
  2102.    va_list arg_marker;  /* variable arg list */
  2103.    int argcount;    /* how many optional args */
  2104.    char *timestring;
  2105.    time_t ltime;
  2106.    FILE *fp;
  2107.    int out;
  2108.    int args[4];
  2109.    int i;
  2110.  
  2111.    /* argc = number of optional arguments */
  2112.    va_start(arg_marker,argc);
  2113.  
  2114.    i=0;
  2115.    while(i < argc)
  2116.       args[i++] = va_arg(arg_marker,int);
  2117.  
  2118.    if(timerflag)
  2119.       fp=fopen("bench","a");
  2120.    timer_start = clock();
  2121.    switch(argc)
  2122.    {
  2123.    case 0:
  2124.       out = (*fractal)();
  2125.       break;
  2126.    case 1:
  2127.       out = (*fractal)(args[0]);
  2128.       break;
  2129.    case 2:
  2130.       out = (*fractal)(args[0],args[1]);
  2131.       break;
  2132.    case 3:
  2133.       out = (*fractal)(args[0],args[1],args[2]);
  2134.       break;
  2135.    case 4:
  2136.       out = (*fractal)(args[0],args[1],args[2],args[3]);
  2137.       break;
  2138.    default:
  2139.       return(-1);
  2140.       break;
  2141.    }
  2142.    /* next assumes CLK_TCK is 10^n, n>=2 */
  2143.    timer_interval = (clock() - timer_start) / (CLK_TCK/100);
  2144.    if(timerflag)
  2145.    {
  2146.       time(<ime);
  2147.       timestring = ctime(<ime);
  2148.       timestring[24] = 0; /*clobber newline in time string */
  2149.  
  2150.       /* at the moment only decode timer takes a parameter */
  2151.       if(argc)
  2152.          fprintf(fp,"decode ");
  2153.       fprintf(fp,"%s type=%s resolution = %dx%d maxiter=%d",
  2154.           timestring,
  2155.           fractalspecific[fractype].name,
  2156.           xdots,
  2157.           ydots,
  2158.           maxit);
  2159.       fprintf(fp," time= %ld.%02ld secs\n",timer_interval/100,timer_interval%100);
  2160.       if(fp != NULL)
  2161.          fclose(fp);
  2162.    }
  2163.    return(out);
  2164. }
  2165.  
  2166. /*
  2167.    I, Timothy Wegner, invented this solidguessing idea and implemented it in
  2168.    more or less the overall framework you see here.  I am adding this note
  2169.    now in a possibly vain attempt to secure my place in history, because
  2170.    Pieter Branderhorst has totally rewritten this routine, incorporating
  2171.    a *MUCH* more sophisticated algorithm.  His revised code is not only
  2172.    faster, but is also more accurate. Harrumph!
  2173. */
  2174.  
  2175. solidguess()
  2176. {
  2177.    int i,x,y,xlim,ylim,blocksize;
  2178.    unsigned int *pfxp0,*pfxp1;
  2179.    unsigned int u;
  2180.  
  2181.    guessplot=(plot!=putcolor && plot!=symplot2 && plot!=symplot2J);
  2182.    /* check if guessing at bottom & right edges is ok */
  2183.    bottom_guess = (plot==symplot2 || (plot==putcolor && iystop+1==ydots));
  2184.    right_guess  = (plot==symplot2J
  2185.        || ((plot==putcolor || plot==symplot2) && ixstop+1==xdots));
  2186.  
  2187.    maxblock=blocksize=ssg_blocksize();
  2188.  
  2189.    calcmode = 1;
  2190.    (*calctype)(); /* initialize calc routine */
  2191.    numpasses = 0; /* for calcmand */
  2192.  
  2193.    /* ensure window top and left are on required boundary, treat window
  2194.          as larger than it really is if necessary (this is the reason symplot
  2195.          routines must check for > xdots/ydots before plotting sym points) */
  2196.    ixstart &= -1 - (maxblock-1);
  2197.    iystart = yybegin;
  2198.    iystart &= -1 - (maxblock-1);
  2199.  
  2200.    if (workpass == 0) /* otherwise first pass already done */
  2201.    {
  2202.       /* first pass, calc every blocksize**2 pixel, quarter result & paint it */
  2203.       if (iystart <= yystart) /* first time for this window, init it */
  2204.       {
  2205.          memset(&prefix[1][0][0],0,maxxblk*maxyblk*2); /* noskip flags off */
  2206.          calcmode=2;
  2207.          row=iystart;
  2208.          for(col=ixstart; col<=ixstop; col+=maxblock)
  2209.          { /* calc top row */
  2210.             if((*calctype)()==-1)
  2211.             {
  2212.                add_worklist(xxstart,xxstop,yystart,yystop,yybegin,0,worksym);
  2213.                goto exit_solidguess;
  2214.             }
  2215.             calcmode = 3;
  2216.          }
  2217.       }
  2218.       else
  2219.          memset(&prefix[1][0][0],-1,maxxblk*maxyblk*2); /* noskip flags on */
  2220.       for(y=iystart; y<=iystop; y+=blocksize)
  2221.       {
  2222.          i = 0;
  2223.          if(y+blocksize<=iystop)
  2224.          { /* calc the row below */
  2225.             row=y+blocksize;
  2226.             calcmode = 2; /* for the first dot on the row */
  2227.             for(col=ixstart; col<=ixstop; col+=maxblock)
  2228.             {
  2229.                if((i=(*calctype)()) == -1)
  2230.                   break;
  2231.                calcmode = 3;
  2232.             }
  2233.          }
  2234.          calcmode=2;
  2235.          if (i == -1 || guessrow(1,y,blocksize) != 0) /* interrupted? */
  2236.          {
  2237.             if (y < yystart)
  2238.                y = yystart;
  2239.             add_worklist(xxstart,xxstop,yystart,yystop,y,0,worksym);
  2240.             goto exit_solidguess;
  2241.          }
  2242.       }
  2243.  
  2244.       if (num_worklist) /* work list not empty, just do 1st pass */
  2245.       {
  2246.          add_worklist(xxstart,xxstop,yystart,yystop,yystart,1,worksym);
  2247.          goto exit_solidguess;
  2248.       }
  2249.       ++workpass;
  2250.       iystart = yystart & (-1 - (maxblock-1));
  2251.  
  2252.       /* calculate skip flags for skippable blocks */
  2253.       xlim=(ixstop+maxblock)/maxblock+1;
  2254.       ylim=((iystop+maxblock)/maxblock+15)/16+1;
  2255.       if(right_guess==0) /* no right edge guessing, zap border */
  2256.          for(y=0;y<=ylim;++y)
  2257.             prefix[1][y][xlim]=-1;
  2258.       if(bottom_guess==0) /* no bottom edge guessing, zap border */
  2259.       {
  2260.          i=(iystop+maxblock)/maxblock+1;
  2261.          y=i/16+1;
  2262.          i=1<<(i&15);
  2263.          for(x=0;x<=xlim;++x)
  2264.             prefix[1][y][x]|=i;
  2265.       }
  2266.       /* set each bit in prefix[0] to OR of it & surrounding 8 in prefix[1] */
  2267.       for(y=0;++y<ylim;)
  2268.       {
  2269.          pfxp0=&prefix[0][y][0];
  2270.          pfxp1=&prefix[1][y][0];
  2271.          for(x=0;++x<xlim;)
  2272.          {
  2273.             ++pfxp1;
  2274.             u=*(pfxp1-1)|*pfxp1|*(pfxp1+1);
  2275.             *(++pfxp0)=u|(u>>1)|(u<<1)
  2276.                |((*(pfxp1-(maxxblk+1))|*(pfxp1-maxxblk)|*(pfxp1-(maxxblk-1)))>>15)
  2277.                   |((*(pfxp1+(maxxblk-1))|*(pfxp1+maxxblk)|*(pfxp1+(maxxblk+1)))<<15);
  2278.          }
  2279.       }
  2280.    }
  2281.    else /* first pass already done */
  2282.       memset(&prefix[0][0][0],-1,maxxblk*maxyblk*2); /* noskip flags on */
  2283.  
  2284.    /* remaining pass(es), halve blocksize & quarter each blocksize**2 */
  2285.    i = workpass;
  2286.    while (--i > 0) /* allow for already done passes */
  2287.       blocksize = blocksize>>1;
  2288.    calcmode=2;
  2289.    while((blocksize=blocksize>>1)>=2)
  2290.    {
  2291.       for(y=iystart; y<=iystop; y+=blocksize)
  2292.          if(guessrow(0,y,blocksize)!=0)
  2293.          {
  2294.             if (y < yystart)
  2295.                y = yystart;
  2296.             add_worklist(xxstart,xxstop,yystart,yystop,y,workpass,worksym);
  2297.             goto exit_solidguess;
  2298.          }
  2299.       ++workpass;
  2300.       if (num_worklist /* work list not empty, do one pass at a time */
  2301.       && blocksize>2) /* if 2, we just did last pass */
  2302.       {
  2303.          add_worklist(xxstart,xxstop,yystart,yystop,yystart,workpass,worksym);
  2304.          goto exit_solidguess;
  2305.       }
  2306.       iystart = yystart & (-1 - (maxblock-1));
  2307.    }
  2308.  
  2309.    exit_solidguess:
  2310.    return(0);
  2311. }
  2312.  
  2313. ssg_blocksize() /* used above and by zoom panning */
  2314. {  
  2315.    int blocksize,i;
  2316.    /* blocksize 4 if <300 rows, 8 if 300-599, 16 if 600-1199, 32 if >=1200 */
  2317.    blocksize=4;
  2318.    i=300;
  2319.    while(i<=ydots)
  2320.    {
  2321.       blocksize+=blocksize;
  2322.       i+=i;
  2323.    }
  2324.    /* increase blocksize if prefix array not big enough */
  2325.    while(blocksize*(maxxblk-2)<xdots || blocksize*(maxyblk-2)*16<ydots)
  2326.       blocksize+=blocksize;
  2327.    return(blocksize);
  2328. }
  2329.  
  2330. #define calcadot(c,x,y) { col=x; row=y; if((c=(*calctype)())==-1) return -1; }
  2331.  
  2332. int guessrow(int firstpass,int y,int blocksize)
  2333. {
  2334.    int x,i,j,color;
  2335.    int xplushalf,xplusblock;
  2336.    int ylessblock,ylesshalf,yplushalf,yplusblock;
  2337.    int     c21,c31,c41;         /* cxy is the color of pixel at (x,y) */
  2338.    int c12,c22,c32,c42;         /* where c22 is the topleft corner of */
  2339.    int c13,c23,c33;             /* the block being handled in current */
  2340.    int     c24,    c44;         /* iteration                          */
  2341.    int guessed23,guessed32,guessed33,guessed12,guessed13;
  2342.    int orig23,orig32,orig33;
  2343.    int prev11,fix21,fix31;
  2344.    unsigned int *pfxptr,pfxmask;
  2345.  
  2346.    halfblock=blocksize>>1;
  2347.    i=y/maxblock;
  2348.    pfxptr=&prefix[firstpass][(i>>4)+1][ixstart/maxblock];
  2349.    pfxmask=1<<(i&15);
  2350.    ylesshalf=y-halfblock;
  2351.    ylessblock=y-blocksize; /* constants, for speed */
  2352.    yplushalf=y+halfblock;
  2353.    yplusblock=y+blocksize;
  2354.    prev11=-1;
  2355.    c24=c12=c13=c22=getcolor(ixstart,y);
  2356.    c31=c21=getcolor(ixstart,(y>0)?ylesshalf:0);
  2357.    if(yplusblock<=iystop)
  2358.       c24=getcolor(ixstart,yplusblock);
  2359.    else if(bottom_guess==0)
  2360.       c24=-1;
  2361.    guessed12=guessed13=0;
  2362.  
  2363.    for(x=ixstart; x<=ixstop;)  /* increment at end, or when doing continue */
  2364.    {
  2365.       if((x&(maxblock-1))==0)  /* time for skip flag stuff */
  2366.       {
  2367.          ++pfxptr;
  2368.          if(firstpass==0 && (*pfxptr&pfxmask)==0)  /* check for fast skip */
  2369.          {
  2370.             /* next useful in testing to make skips visible */
  2371.             /*
  2372.                            if(halfblock==1)
  2373.                            {
  2374.                               (*plot)(x+1,y,0); (*plot)(x,y+1,0); (*plot)(x+1,y+1,0);
  2375.                               }
  2376.                          */
  2377.             x+=maxblock;
  2378.             prev11=c31=c21=c24=c12=c13=c22;
  2379.             guessed12=guessed13=0;
  2380.             continue;
  2381.          }
  2382.       }
  2383.  
  2384.       if(firstpass)  /* 1st pass, paint topleft corner */
  2385.          plotblock(0,x,y,c22);
  2386.       /* setup variables */
  2387.       xplusblock=(xplushalf=x+halfblock)+halfblock;
  2388.       if(xplushalf>ixstop)
  2389.       {
  2390.          if(right_guess==0)
  2391.             c31=-1;
  2392.       }
  2393.       else if(y>0)
  2394.          c31=getcolor(xplushalf,ylesshalf);
  2395.       if(xplusblock<=ixstop)
  2396.       {
  2397.          if(yplusblock<=iystop)
  2398.             c44=getcolor(xplusblock,yplusblock);
  2399.          c41=getcolor(xplusblock,(y>0)?ylesshalf:0);
  2400.          c42=getcolor(xplusblock,y);
  2401.       }
  2402.       else if(right_guess==0)
  2403.          c41=c42=c44=-1;
  2404.       if(yplusblock>iystop)
  2405.          c44=(bottom_guess)?c42:-1;
  2406.  
  2407.       /* guess or calc the remaining 3 quarters of current block */
  2408.       guessed23=guessed32=guessed33=1;
  2409.       c23=c32=c33=c22;
  2410.       if(yplushalf>iystop)
  2411.       {
  2412.          if(bottom_guess==0)
  2413.             c23=c33=-1;
  2414.          guessed23=guessed33=-1;
  2415.       }
  2416.       if(xplushalf>ixstop)
  2417.       {
  2418.          if(right_guess==0)
  2419.             c32=c33=-1;
  2420.          guessed32=guessed33=-1;
  2421.       }
  2422.       while(1) /* go around till none of 23,32,33 change anymore */
  2423.       {
  2424.          if(guessed33>0
  2425.              && (c33!=c44 || c33!=c42 || c33!=c24 || c33!=c32 || c33!=c23))
  2426.          {
  2427.             calcadot(c33,xplushalf,yplushalf);
  2428.             guessed33=0;
  2429.          }
  2430.          if(guessed32>0
  2431.              && (c32!=c33 || c32!=c42 || c32!=c31 || c32!=c21
  2432.              || c32!=c41 || c32!=c23))
  2433.          {
  2434.             calcadot(c32,xplushalf,y);
  2435.             guessed32=0;
  2436.             continue;
  2437.          }
  2438.          if(guessed23>0
  2439.              && (c23!=c33 || c23!=c24 || c23!=c13 || c23!=c12 || c23!=c32))
  2440.          {
  2441.             calcadot(c23,x,yplushalf);
  2442.             guessed23=0;
  2443.             continue;
  2444.          }
  2445.          break;
  2446.       }
  2447.  
  2448.       if(firstpass) /* note whether any of block's contents were calculated */
  2449.          if(guessed23==0 || guessed32==0 || guessed33==0)
  2450.             *pfxptr|=pfxmask;
  2451.  
  2452.       if(halfblock>1) /* not last pass, check if something to display */
  2453.          if(firstpass)  /* display guessed corners, fill in block */
  2454.          {
  2455.             if(guessplot)
  2456.             {
  2457.                if(guessed23>0)
  2458.                   (*plot)(x,yplushalf,c23);
  2459.                if(guessed32>0)
  2460.                   (*plot)(xplushalf,y,c32);
  2461.                if(guessed33>0)
  2462.                   (*plot)(xplushalf,yplushalf,c33);
  2463.             }
  2464.             plotblock(1,x,yplushalf,c23);
  2465.             plotblock(0,xplushalf,y,c32);
  2466.             plotblock(1,xplushalf,yplushalf,c33);
  2467.          }
  2468.          else  /* repaint changed blocks */
  2469.          {
  2470.             if(c23!=c22)
  2471.                plotblock(-1,x,yplushalf,c23);
  2472.             if(c32!=c22)
  2473.                plotblock(-1,xplushalf,y,c32);
  2474.             if(c33!=c22)
  2475.                plotblock(-1,xplushalf,yplushalf,c33);
  2476.          }
  2477.  
  2478.       /* check if some calcs in this block mean earlier guesses need fixing */
  2479.       fix21=((c22!=c12 || c22!=c32)
  2480.           && c21==c22 && c21==c31 && c21==prev11
  2481.           && y>0
  2482.           && (x==ixstart || c21==getcolor(x-halfblock,ylessblock))
  2483.           && (xplushalf>ixstop || c21==getcolor(xplushalf,ylessblock))
  2484.           && c21==getcolor(x,ylessblock));
  2485.       fix31=(c22!=c32
  2486.           && c31==c22 && c31==c42 && c31==c21 && c31==c41
  2487.           && y>0 && xplushalf<=ixstop
  2488.           && c31==getcolor(xplushalf,ylessblock)
  2489.           && (xplusblock>ixstop || c31==getcolor(xplusblock,ylessblock))
  2490.           && c31==getcolor(x,ylessblock));
  2491.       prev11=c31; /* for next time around */
  2492.       if(fix21)
  2493.       {
  2494.          calcadot(c21,x,ylesshalf);
  2495.          if(halfblock>1 && c21!=c22)
  2496.             plotblock(-1,x,ylesshalf,c21);
  2497.       }
  2498.       if(fix31)
  2499.       {
  2500.          calcadot(c31,xplushalf,ylesshalf);
  2501.          if(halfblock>1 && c31!=c22)
  2502.             plotblock(-1,xplushalf,ylesshalf,c31);
  2503.       }
  2504.       if(c23!=c22)
  2505.       {
  2506.          if(guessed12)
  2507.          {
  2508.             calcadot(c12,x-halfblock,y);
  2509.             if(halfblock>1 && c12!=c22)
  2510.                plotblock(-1,x-halfblock,y,c12);
  2511.          }
  2512.          if(guessed13)
  2513.          {
  2514.             calcadot(c13,x-halfblock,yplushalf);
  2515.             if(halfblock>1 && c13!=c22)
  2516.                plotblock(-1,x-halfblock,yplushalf,c13);
  2517.          }
  2518.       }
  2519.       c22=c42;
  2520.       c24=c44;
  2521.       c13=c33;
  2522.       c31=c21=c41;
  2523.       c12=c32;
  2524.       guessed12=guessed32;
  2525.       guessed13=guessed33;
  2526.       x+=blocksize;
  2527.    } /* end x loop */
  2528.  
  2529.    if(firstpass==0 || guessplot) return 0;
  2530.  
  2531.    /* paint rows the fast way */
  2532.    for(i=0;i<halfblock;++i)
  2533.    {
  2534.       if((j=y+i)<=iystop)
  2535.          put_line(j,xxstart,ixstop,&dstack[xxstart]);
  2536.       if((j=y+i+halfblock)<=iystop)
  2537.          put_line(j,xxstart,ixstop,&dstack[xxstart+2048]);
  2538.    }
  2539.    if(plot!=putcolor)  /* symmetry, just vertical & origin the fast way */
  2540.    {
  2541.       if(plot==symplot2J) /* origin sym, reverse lines */
  2542.          for(i=(ixstop+xxstart+1)/2;--i>=xxstart;)
  2543.          {
  2544.             color=dstack[i];
  2545.             dstack[i]=dstack[j=ixstop-(i-xxstart)];
  2546.             dstack[j]=color;
  2547.             j+=2048;
  2548.             color=dstack[i+2048];
  2549.             dstack[i+2048]=dstack[j];
  2550.             dstack[j]=color;
  2551.          }
  2552.       for(i=0;i<halfblock;++i)
  2553.       {
  2554.          if((j=yystop-(y+i-yystart))>iystop && j<ydots)
  2555.             put_line(j,xxstart,ixstop,&dstack[xxstart]);
  2556.          if((j=yystop-(y+i+halfblock-yystart))>iystop && j<ydots)
  2557.             put_line(j,xxstart,ixstop,&dstack[xxstart+2048]);
  2558.       }
  2559.    }
  2560.    return 0;
  2561. }
  2562.  
  2563. plotblock(int buildrow,int x,int y,int color)
  2564. {
  2565.    int i,xlim,ylim;
  2566.    if((xlim=x+halfblock)>ixstop)
  2567.       xlim=ixstop+1;
  2568.    if(buildrow>=0 && guessplot==0) /* save it for later put_line */
  2569.    {
  2570.       if(buildrow==0)
  2571.          for(i=x;i<xlim;++i)
  2572.             dstack[i]=color;
  2573.       else
  2574.          for(i=x;i<xlim;++i)
  2575.             dstack[i+2048]=color;
  2576.       if (x>=xxstart) /* when x reduced for alignment, paint those dots too */
  2577.          return(0); /* the usual case */
  2578.    }
  2579.    /* paint it */
  2580.    if((ylim=y+halfblock)>iystop)
  2581.    {
  2582.       if(y>iystop)
  2583.          return(0);
  2584.       ylim=iystop+1;
  2585.    }
  2586.    for(i=x;++i<xlim;)
  2587.       (*plot)(i,y,color); /* skip 1st dot on 1st row */
  2588.    while(++y<ylim)
  2589.       for(i=x;i<xlim;++i)
  2590.          (*plot)(i,y,color);
  2591. }
  2592.  
  2593. /******************************************************************/
  2594. /* Continuous potential calculation for Mandelbrot and Julia      */
  2595. /* Reference: Science of Fractal Images p. 190.                   */
  2596. /* Special thanks to Mark Peterson for his "MtMand" program that  */
  2597. /* beautifully approximates plate 25 (same reference) and spurred */
  2598. /* on the inclusion of similar capabilities in FRACTINT.          */
  2599. /*                                                                */
  2600. /* The purpose of this function is to calculate a color value     */
  2601. /* for a fractal that varies continuously with the screen pixels  */
  2602. /* locations for better rendering in 3D.                          */
  2603. /*                                                                */
  2604. /* Here "magnitude" is the modulus of the orbit value at          */
  2605. /* "iterations". The potparms[] are user-entered paramters        */
  2606. /* controlling the level and slope of the continuous potential    */
  2607. /* surface. Returns color.  - Tim Wegner 6/25/89                  */
  2608. /*                                                                */
  2609. /*                     -- Change history --                       */
  2610. /*                                                                */
  2611. /* 09/12/89   - added floatflag support and fixed float underflow */
  2612. /*                                                                */
  2613. /******************************************************************/
  2614.  
  2615. int potential(double mag, int iterations)
  2616. {
  2617.    extern char potfile[];     /* name of pot save file */
  2618.    extern struct fractal_info save_info;        /*  for saving data */
  2619.  
  2620.    static int x,y;            /* keep track of position in image */
  2621.    float f_mag,f_tmp,pot;
  2622.    double d_tmp;
  2623.    int i_pot;
  2624.    unsigned int intbuf;
  2625.    FILE *t16_create();
  2626.  
  2627.    extern int maxit;          /* maximum iteration bailout limit */
  2628.    extern double potparam[];     /* continuous potential parameters */
  2629.    extern char potfile[];     /* pot savename */
  2630.    extern char floatflag;
  2631.  
  2632.    if(top)
  2633.    {
  2634.       top = 0;
  2635.       x   = 0;
  2636.       y   = 0;
  2637.       /* make sure file is closed */
  2638.       if(fp_pot)
  2639.          fclose(fp_pot);
  2640.       if(potfile[0])
  2641.          fp_pot = t16_create(potfile, xdots, ydots,
  2642.              sizeof(struct fractal_info), &save_info);
  2643.    }
  2644.    if(iterations < maxit)
  2645.    {
  2646.       i_pot = pot = iterations+2;
  2647.       if(i_pot <= 0 || mag <= 1.0)
  2648.          pot = 0.0;
  2649.       else /* pot = log(mag) / pow(2.0, (double)pot); */
  2650.       {
  2651.          i_pot = pot;
  2652.          if(i_pot < 120 && !floatflag) /* empirically determined limit of fShift */
  2653.          {
  2654.             f_mag = mag;
  2655.             fLog14(f_mag,f_tmp); /* this SHOULD be non-negative */
  2656.             fShift(f_tmp,-i_pot,pot);
  2657.          }
  2658.          else
  2659.          {
  2660.             d_tmp = log(mag)/(double)pow(2.0,(double)pot);
  2661.             if(d_tmp > FLT_MIN) /* prevent float type underflow */
  2662.                pot = d_tmp;
  2663.             else
  2664.                pot = 0.0;
  2665.          }
  2666.       }
  2667.       /* following transformation strictly for aesthetic reasons */
  2668.       /* meaning of parameters:
  2669.             potparam[0] -- zero potential level - highest color -
  2670.             potparam[1] -- slope multiplier -- higher is steeper
  2671.             potparam[2] -- rqlim value if changeable (bailout for modulus) */
  2672.  
  2673.       if(pot > 0.0)
  2674.       {
  2675.          if(floatflag)
  2676.             pot = (float)sqrt((double)pot);
  2677.          else
  2678.          {
  2679.             fSqrt14(pot,f_tmp);
  2680.             pot = f_tmp;
  2681.          }
  2682.          pot = potparam[0] - pot*potparam[1] - 1.0;
  2683.       }
  2684.       else
  2685.          pot = potparam[0] - 1.0;
  2686.    }
  2687.    else if(inside > 0)
  2688.       pot = inside; /* negative inside turns off inside replacement */
  2689.    else
  2690.       pot = potparam[0];
  2691.  
  2692.    if(pot > colors)
  2693.       pot = colors -1;
  2694.    if(pot < 1.0)
  2695.       pot = 1.0; /* avoid color 0 */
  2696.  
  2697.    if(fp_pot)
  2698.    {
  2699.       intbuf = pot*(1<<8); /* shift 8 bits */
  2700.       boxx[x] = intbuf;
  2701.       if(++x == xdots)
  2702.       {
  2703.          /* new row */
  2704.          t16_putline(fp_pot, xdots, boxx);
  2705.          x = 0;
  2706.          if(++y == ydots)
  2707.          {
  2708.             /* we're done */
  2709.             t16_flush(fp_pot);
  2710.             fclose(fp_pot);
  2711.             fp_pot = NULL;
  2712.             potfile[0] = 0;
  2713.          }
  2714.       }
  2715.    }
  2716.    return((int)pot);
  2717. }
  2718.  
  2719. /*
  2720.    "intpotential" is called by the MANDEL/JULIA routines with scaled long
  2721.    integer magnitudes.  The routine just calls "potential".  Bert
  2722. */
  2723.  
  2724. int intpotential(unsigned long longmagnitude, int iterations)
  2725. {                                /* this version called from calcmand() */
  2726.    double magnititude;
  2727.  
  2728.    magnitude = ((double)longmagnitude) / fudge;
  2729.    magnitude = 256*magnitude;
  2730.    return(potential(magnitude, iterations));
  2731. }
  2732.  
  2733. static int xsym_split(int xaxis_row,int xaxis_between)
  2734. {
  2735.    int i;
  2736.    if ((worksym&0x11) == 0x10) /* already decided not sym */
  2737.       return(1);
  2738.    if ((worksym&1) != 0) /* already decided on sym */
  2739.       iystop = (yystart+yystop)/2;
  2740.    else /* new window, decide */
  2741.    {
  2742.       worksym |= 0x10;
  2743.       if (xaxis_row <= yystart || xaxis_row >= yystop)
  2744.          return(1); /* axis not in window */
  2745.       i = xaxis_row + (xaxis_row - yystart);
  2746.       if (xaxis_between)
  2747.          ++i;
  2748.       if (i > yystop) /* split into 2 pieces, bottom has the symmetry */
  2749.       {
  2750.          if (num_worklist >= MAXCALCWORK-1) /* no room to split */
  2751.             return(1);
  2752.          iystop = xaxis_row - (yystop - xaxis_row);
  2753.          if (!xaxis_between)
  2754.             --iystop;
  2755.          add_worklist(xxstart,xxstop,iystop+1,yystop,iystop+1,workpass,0);
  2756.          yystop = iystop;
  2757.          return(1); /* tell set_symmetry no sym for current window */
  2758.       }
  2759.       if (i < yystop) /* split into 2 pieces, top has the symmetry */
  2760.       {
  2761.          if (num_worklist >= MAXCALCWORK-1) /* no room to split */
  2762.             return(1);
  2763.          add_worklist(xxstart,xxstop,i+1,yystop,i+1,workpass,0);
  2764.          yystop = i;
  2765.       }
  2766.       iystop = xaxis_row;
  2767.       worksym |= 1;
  2768.    }
  2769.    symmetry = 0;
  2770.    return(0); /* tell set_symmetry its a go */
  2771. }
  2772.  
  2773. static int ysym_split(int yaxis_col,int yaxis_between)
  2774. {
  2775.    int i;
  2776.    if ((worksym&0x22) == 0x20) /* already decided not sym */
  2777.       return(1);
  2778.    if ((worksym&2) != 0) /* already decided on sym */
  2779.       ixstop = (xxstart+xxstop)/2;
  2780.    else /* new window, decide */
  2781.    {
  2782.       worksym |= 0x20;
  2783.       if (yaxis_col <= xxstart || yaxis_col >= xxstop)
  2784.          return(1); /* axis not in window */
  2785.       i = yaxis_col + (yaxis_col - xxstart);
  2786.       if (yaxis_between)
  2787.          ++i;
  2788.       if (i > xxstop) /* split into 2 pieces, right has the symmetry */
  2789.       {
  2790.          if (num_worklist >= MAXCALCWORK-1) /* no room to split */
  2791.             return(1);
  2792.          ixstop = yaxis_col - (xxstop - yaxis_col);
  2793.          if (!yaxis_between)
  2794.             --ixstop;
  2795.          add_worklist(ixstop+1,xxstop,yystart,yystop,yystart,workpass,0);
  2796.          xxstop = ixstop;
  2797.          return(1); /* tell set_symmetry no sym for current window */
  2798.       }
  2799.       if (i < xxstop) /* split into 2 pieces, left has the symmetry */
  2800.       {
  2801.          if (num_worklist >= MAXCALCWORK-1) /* no room to split */
  2802.             return(1);
  2803.          add_worklist(i+1,xxstop,yystart,yystop,yystart,workpass,0);
  2804.          xxstop = i;
  2805.       }
  2806.       ixstop = yaxis_col;
  2807.       worksym |= 2;
  2808.    }
  2809.    symmetry = 0;
  2810.    return(0); /* tell set_symmetry its a go */
  2811. }
  2812.  
  2813. static void setsymmetry(int sym, int uselist) /* set up proper symmetrical plot functions */
  2814. {
  2815.    extern int forcesymmetry;
  2816.    int parmszero;
  2817.    int xaxis_row, yaxis_col;         /* pixel number for origin */
  2818.    int xaxis_between, yaxis_between; /* if axis between 2 pixels, not on one */
  2819.    double ftemp;
  2820.    symmetry = 1;
  2821.    if(sym == NOPLOT && forcesymmetry == 999)
  2822.    {
  2823.       plot = noplot;
  2824.       return;
  2825.    }
  2826.    /* NOTE: 16-bit potential disables symmetry */
  2827.    /* also any decomp= option and any inversion not about the origin */
  2828.    /* also any rotation other than 180deg and any off-axis stretch */
  2829.    if (potfile[0] || (invert && inversion[2] != 0.0)
  2830.        || decomp[0] != 0
  2831.        || xxmin!=xx3rd || yymin!=yy3rd)
  2832.       return;
  2833.    if(sym != XAXIS && sym != XAXIS_NOPARM && inversion[1] != 0.0 && forcesymmetry == 999)
  2834.       return;
  2835.    if(forcesymmetry != 999)
  2836.       sym = forcesymmetry;
  2837.    parmszero = (parm.x == 0.0 && parm.y == 0.0 && useinitorbit != 1);
  2838.    xaxis_row = yaxis_col = -1;
  2839.    if (fabs(yymin+yymax) < fabs(yymin)+fabs(yymax)) /* axis is on screen */
  2840.    {
  2841.       ftemp = (0.0-yymax) / (yymin-yymax) * (ydots-1) + 0.25;
  2842.       xaxis_row = ftemp;
  2843.       xaxis_between = (ftemp - xaxis_row >= 0.5);
  2844.       if (uselist == 0 && (!xaxis_between || (xaxis_row+1)*2 != ydots))
  2845.          xaxis_row = -1; /* can't split screen, so dead center or not at all */
  2846.    }
  2847.    if (fabs(xxmin+xxmax) < fabs(xxmin)+fabs(xxmax)) /* axis is on screen */
  2848.    {
  2849.       ftemp = (0.0-xxmin) / (xxmax-xxmin) * (xdots-1) + 0.25;
  2850.       yaxis_col = ftemp;
  2851.       yaxis_between = (ftemp - yaxis_col >= 0.5);
  2852.       if (uselist == 0 && (!yaxis_between || (yaxis_col+1)*2 != xdots))
  2853.          yaxis_col = -1; /* can't split screen, so dead center or not at all */
  2854.    }
  2855.    switch(sym)       /* symmetry switch */
  2856.    {
  2857.    case XAXIS_NOREAL:    /* X-axis Symmetry (no real param) */
  2858.       if (parm.x != 0.0) break;
  2859.       goto xsym;
  2860.    case XAXIS_NOIMAG:    /* X-axis Symmetry (no imag param) */
  2861.       if (parm.y != 0.0) break;
  2862.       goto xsym;
  2863.    case XAXIS_NOPARM:                        /* X-axis Symmetry  (no params)*/
  2864.       if (!parmszero)
  2865.          break;
  2866.       xsym:
  2867.    case XAXIS:                       /* X-axis Symmetry */
  2868.       if (xsym_split(xaxis_row,xaxis_between) == 0)
  2869.          if(basin)
  2870.             plot = symplot2basin;
  2871.          else
  2872.             plot = symplot2;
  2873.       break;
  2874.    case YAXIS_NOPARM:                        /* Y-axis Symmetry (No Parms)*/
  2875.       if (!parmszero)
  2876.          break;
  2877.    case YAXIS:                       /* Y-axis Symmetry */
  2878.       if (ysym_split(yaxis_col,yaxis_between) == 0)
  2879.          plot = symplot2Y;
  2880.       break;
  2881.    case XYAXIS_NOPARM:                       /* X-axis AND Y-axis Symmetry (no parms)*/
  2882.       if(!parmszero)
  2883.          break;
  2884.    case XYAXIS:                      /* X-axis AND Y-axis Symmetry */
  2885.       xsym_split(xaxis_row,xaxis_between);
  2886.       ysym_split(yaxis_col,yaxis_between);
  2887.       switch (worksym & 3)
  2888.       {
  2889.       case 1: /* just xaxis symmetry */
  2890.          if(basin)
  2891.             plot = symplot2basin;
  2892.          else
  2893.             plot = symplot2 ;
  2894.          break;
  2895.       case 2: /* just yaxis symmetry */
  2896.          if (basin) /* got no routine for this case */
  2897.          {
  2898.             ixstop = xxstop; /* fix what split should not have done */
  2899.             symmetry = 1;
  2900.          }
  2901.          else
  2902.             plot = symplot2Y;
  2903.          break;
  2904.       case 3: /* both axes */
  2905.          if(basin)
  2906.             plot = symplot4basin;
  2907.          else
  2908.             plot = symplot4 ;
  2909.       }
  2910.       break;
  2911.    case ORIGIN_NOPARM:                       /* Origin Symmetry (no parms)*/
  2912.       if (!parmszero)
  2913.          break;
  2914.    case ORIGIN:                      /* Origin Symmetry */
  2915.       originsym:
  2916.       if ( xsym_split(xaxis_row,xaxis_between) == 0
  2917.           && ysym_split(yaxis_col,yaxis_between) == 0)
  2918.       {
  2919.          plot = symplot2J;
  2920.          ixstop = xxstop; /* didn't want this changed */
  2921.       }
  2922.       else
  2923.       {
  2924.          iystop = yystop; /* in case first split worked */
  2925.          symmetry = 1;
  2926.          worksym = 0x30; /* let it recombine with others like it */
  2927.       }
  2928.       break;
  2929.    case PI_SYM_NOPARM:
  2930.       if (!parmszero)
  2931.          break;
  2932.    case PI_SYM:                      /* PI symmetry */
  2933.       if(invert && forcesymmetry == 999)
  2934.          goto originsym;
  2935.       plot = symPIplot ;
  2936.       symmetry = 0;
  2937.       if ( xsym_split(xaxis_row,xaxis_between) == 0
  2938.           && ysym_split(yaxis_col,yaxis_between) == 0)
  2939.          if(parm.y == 0.0)
  2940.             plot = symPIplot4J; /* both axes */
  2941.          else
  2942.             plot = symPIplot2J; /* origin */
  2943.       else
  2944.       {
  2945.          iystop = yystop; /* in case first split worked */
  2946.          worksym = 0x30;  /* don't mark pisym as ysym, just do it unmarked */
  2947.       }
  2948.       pixelpi = (PI/fabs(xxmax-xxmin))*xdots; /* PI in pixels */
  2949.       if ( (ixstop = xxstart+pixelpi-1 ) > xxstop)
  2950.          ixstop = xxstop;
  2951.       if (plot == symPIplot4J && ixstop > (i = (xxstart+xxstop)/2))
  2952.          ixstop = i;
  2953.       break;
  2954.    default:                  /* no symmetry */
  2955.       break;
  2956.    }
  2957. }
  2958.  
  2959. int add_worklist(int xfrom, int xto, int yfrom, int yto, int ybegin,
  2960. int pass, int sym)
  2961. {
  2962.    int i;
  2963.    if (num_worklist >= MAXCALCWORK)
  2964.       return(-1);
  2965.    worklist[num_worklist].xxstart = xfrom;
  2966.    worklist[num_worklist].xxstop  = xto;
  2967.    worklist[num_worklist].yystart = yfrom;
  2968.    worklist[num_worklist].yystop  = yto;
  2969.    worklist[num_worklist].yybegin = ybegin;
  2970.    worklist[num_worklist].pass    = pass;
  2971.    worklist[num_worklist].sym     = sym;
  2972.    ++num_worklist;
  2973.    tidy_worklist();
  2974.    return(0);
  2975. }
  2976.  
  2977. static int combine_worklist() /* look for 2 entries which can freely merge */
  2978. {
  2979.    int i,j;
  2980.    for (i=0; i<num_worklist; ++i)
  2981.       if (worklist[i].yystart == worklist[i].yybegin)
  2982.          for (j=i+1; j<num_worklist; ++j)
  2983.             if (worklist[j].sym == worklist[i].sym
  2984.                 && worklist[j].yystart == worklist[j].yybegin
  2985.                 && worklist[i].pass == worklist[j].pass)
  2986.             {
  2987.                if ( worklist[i].xxstart == worklist[j].xxstart
  2988.                    && worklist[i].xxstop  == worklist[j].xxstop)
  2989.                {
  2990.                   if (worklist[i].yystop+1 == worklist[j].yystart)
  2991.                   {
  2992.                      worklist[i].yystop = worklist[j].yystop;
  2993.                      return(j);
  2994.                   }
  2995.                   if (worklist[j].yystop+1 == worklist[i].yystart)
  2996.                   {
  2997.                      worklist[i].yystart = worklist[j].yystart;
  2998.                      worklist[i].yybegin = worklist[j].yybegin;
  2999.                      return(j);
  3000.                   }
  3001.                }
  3002.                if ( worklist[i].yystart == worklist[j].yystart
  3003.                    && worklist[i].yystop  == worklist[j].yystop)
  3004.                {
  3005.                   if (worklist[i].xxstop+1 == worklist[j].xxstart)
  3006.                   {
  3007.                      worklist[i].xxstop = worklist[j].xxstop;
  3008.                      return(j);
  3009.                   }
  3010.                   if (worklist[j].xxstop+1 == worklist[i].xxstart)
  3011.                   {
  3012.                      worklist[i].xxstart = worklist[j].xxstart;
  3013.                      return(j);
  3014.                   }
  3015.                }
  3016.             }
  3017.    return(0); /* nothing combined */
  3018. }
  3019.  
  3020. void tidy_worklist() /* combine mergeable entries, resort */
  3021. {
  3022.    int i,j;
  3023.    struct workliststuff tempwork;
  3024.    while (i=combine_worklist())
  3025.    { /* merged two, delete the gone one */
  3026.       while (++i < num_worklist)
  3027.          worklist[i-1] = worklist[i];
  3028.       --num_worklist;
  3029.    }
  3030.    for (i=0; i<num_worklist; ++i)
  3031.       for (j=i+1; j<num_worklist; ++j)
  3032.          if (worklist[j].pass < worklist[i].pass
  3033.              || (worklist[j].pass == worklist[i].pass
  3034.              && (worklist[j].yystart < worklist[i].yystart
  3035.              || (   worklist[j].yystart == worklist[i].yystart
  3036.              && worklist[j].xxstart <  worklist[i].xxstart))))
  3037.          { /* dumb sort, swap 2 entries to correct order */
  3038.             tempwork = worklist[i];
  3039.             worklist[i] = worklist[j];
  3040.             worklist[j] = tempwork;
  3041.          }
  3042. }
  3043.  
  3044. void get_julia_attractor (double real, double imag)
  3045. {
  3046.    struct lcomplex lresult;
  3047.    struct complex result;
  3048.    int savper,savmaxit;
  3049.  
  3050.    if (attractors == 0 && finattract == 0) /* not magnet & not requested */
  3051.       return;
  3052.  
  3053.    if (attractors >= N_ATTR)     /* space for more attractors ?  */
  3054.       return;                  /* Bad luck - no room left !    */
  3055.  
  3056.    savper = periodicitycheck;
  3057.    savmaxit = maxit;
  3058.    periodicitycheck = 0;
  3059.    old.x = real;                    /* prepare for f.p orbit calc */
  3060.    old.y = imag;
  3061.    tempsqrx = sqr(old.x);
  3062.    tempsqry = sqr(old.y);
  3063.  
  3064.    lold.x = real;           /* prepare for int orbit calc */
  3065.    lold.y = imag;
  3066.    ltempsqrx = tempsqrx;
  3067.    ltempsqry = tempsqry;
  3068.  
  3069.    lold.x = lold.x << bitshift;
  3070.    lold.y = lold.y << bitshift;
  3071.    ltempsqrx = ltempsqrx << bitshift;
  3072.    ltempsqry = ltempsqry << bitshift;
  3073.  
  3074.    if (maxit < 500)         /* we're going to try at least this hard */
  3075.       maxit = 500;
  3076.    color = 0;
  3077.    while (++color < maxit)
  3078.       if(fractalspecific[fractype].orbitcalc())
  3079.          break;
  3080.    if (color >= maxit)      /* if orbit stays in the lake */
  3081.    {
  3082.       if (integerfractal)   /* remember where it went to */
  3083.          lresult = lnew;
  3084.       else
  3085.          result =  new;
  3086.       if(!fractalspecific[fractype].orbitcalc()) /* if it stays in the lake */
  3087.       {                        /* and doen't move far, probably */
  3088.          if (integerfractal)   /*   found a finite attractor    */
  3089.          {
  3090.             if(labs(lresult.x-lnew.x) < lclosenuff
  3091.                 && labs(lresult.y-lnew.y) < lclosenuff)
  3092.             {
  3093.                lattr[attractors] = lnew;
  3094.                attractors++;   /* another attractor - coloured lakes ! */
  3095.             }
  3096.          }
  3097.          else
  3098.          {
  3099.             if(fabs(result.x-new.x) < closenuff
  3100.                 && fabs(result.y-new.y) < closenuff)
  3101.             {
  3102.                attr[attractors] = new;
  3103.                attractors++;   /* another attractor - coloured lakes ! */
  3104.             }
  3105.          }
  3106.       }
  3107.    }
  3108.    if(attractors==0)
  3109.       periodicitycheck = savper;
  3110.    maxit = savmaxit;
  3111. }
  3112.  
  3113.